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PROGRAM USER’S MANUAL FOR OPTIMIZING THE DESIGN OF A 
LIQUID OR GASEOUS PROPELLANT ROCKET ENGINE WITH THE 
AUTOMATED COMBUSTOR DESIGN CODE AUTOCOM 

R. H. Reichel, D. S. Hague, R. T. Jones, and C. R. Glatt 
PART I - DATA INPUT AND OUTPUT WITH A SAMPLE ENGINE TEST CASE 

. >r. ■ • ■ 

r 1. INTRODUCTION 

In designing a liquid or gaseous propellant rocket combustion chamber the 
engineer must compromise between characteristics such as performance, stability, 
weight, injector complexity, cost, etc. These engine aharaateristids are not 
items which are directly controlled by the designer. Instead, they are compli- 
cated functions of the independent design variables available to the designer, 
for example, injector hole size, chamber length, etc. To further complicate the 
problem, frequently there are several techniques that can be used to predict how 
an engine characteristic (such as performance) varies with the independent 
design variables. An engine characteristic value computed by any one of these 
techniques will be referred to as a specific engine characteristic. 

If the engine designer had infinite funds and time available to him, he 
could design many combustors with different combinations and permutations of 
the various independent design variables. Specific engine characteristics could 
then be calculated for each design with all the available techniques if the 
designer had the ability to digest his particular application.. The selection 
would be made, of necessity, on the basis of weighting factors. These weighting 
factors would first be used to define an average combustor characteristic for 
each significant engine characteristic such as performance, stability, etc. 

These average combustor characteristics would then be weighted into a single 
combustor value function or rating which could be used, for example, as a measure 
of the payload loss associated with an ideal engine. 

With limited funds and time, the designer can only examine , a few designs, 
and, because he is familiar with only a few techniques for calculating the 
various characteristics, he uses only this limited set of techniques to test 
the acceptability of each design. Using this approach some characteristics 
are never determined until after the combustor has been built, tested, and often 
found unacceptable. For example, stability characteristics which are partic- 
ularly difficult to assess frequently result in an unacceptable engine design. 
Usually the designs selected in a project are those that are very similar to 
designs that have been successful in the past. As a result, a design of another 
group that would be better for a particular application is frequently neglected 
or ignored. Similarly, when trouble is encountered with a particular engine 
characteristic during the development phase, changes are made to overcome the 
particular problem using past experience. An improved engine design probably 
could have been obtained by determining which design variable or set of vari- 
ables should be used to overcome the problem with the least sacrifice to other 
characteristics . 
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The work performed under the present contract was directed to the 
development of a generalized computer program to calculate dll the character-, 
istics of a given , combustor design. The program then uses a perturbation . 
technique to determine the changes in the design variables that produce the 
greatest improvement in the rating of the combustoi^ design. The program then 
follows sequentially the design path that produces the greatest improvement 
in the rating to arrive at a combustor design that has the best combination 
of all variables. This design as called the optimum combustor design. The ,, 
automated combustor design code which generates the optimum combustor design^ 
has been given the acronym AUTOCOM. An interim report describing the design 
code AUTOCOM was issued earlier. Reference 1-1. 

In any optimization situation, the engineer-designer is faced with the 
problem of selecting the. rating or value function which is to be minimized. 

In this report the rating of a design is based on a weighted average of ^11 
the characteristics of a given combustor. Weighting factors are used to obtain 
both average characteristics and a unique combustor rating. The constants 
first allow the designer to introduce his views regarding the importance or 
validity of one technique for obtaining a given characteristic versus another 
technique for obtaining the same characteristic. For example, if the designer 
believes that only one technique is valid for predicting the performance of 
a given design, he will assign a unity value to the weighting factor constants 
for that specific characteristic, and the constants for all the other perfor- 
mance characteristics will be zero. Second, the weighting factor constants 
in the equation to obtain a single rating for a given combustor design allow 
the designer to introduce the relative importance of different types of 
characteristics, for example, stability versus performance. The weighting 
factor constants, therefore, give the designer the same control and flexibility 
in the corrputer program as he has in the present cut- and- try** system. To 
establish a base point for the rating system, a hypothetical ideal combustor 
is given a rating of zero. The hypothetical ideal combustor would, have one 
hundred per cent of theoretical C* performance, infinite damping rate for all 
modes of instability, zero pressure drop, zero chamber length, chamber diameter 
equal to the throat diameter, etc. 

Techniques for obtaining the various specific engine and average combustor 
characteristics contained in the AUTOCOM code are outlined in this report. 

The code is written in a modular fashion which permits rapid extension of the 
combustor characteristic equations. This approach leads to an open ended 
engine model capable of future development and extension consistent with the 
growth of capability in combustor design analysis. 

Designs are optimized by nonlinear multivariable search techniques embodied 
in an optimization program module. A variety of multivariable search techniques 
may be employed to define an optimal design including 

a. One-parameter-at-a-time techniques 

b. Organized techniques based on the computation of first- or 
second-order rating derivatives with respect to the design 
variables 

c. Randomized techniques 
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Program operation is shown conceptually in Figure 1-1. Fixed engine and 
propellant data is input on punched cards together with the free design 
variable values. The specific engine characteristics are then computed and 
combined through user-specified weighting constants to form the average 
combustor characteristics for performance, stability, etc. The average 
combustor characteristics are then combined through a second set of user- 
specified weighting factors into the combustor rating which measures payload 
loss based on an ideal engine. The rated engine free design parameters are 
then perturbed in the optimization loop until an optimal design is evolved. ‘ 

For the user's convenience this document is organized in two parts. 

Part I deals with the input data setup, the program outputs, and concludes 
with an engine test case. Part II represents the program backup and discusses 
the combustor program structure, the programming techniques, and the AUTOCOM 
program ailalysis. ' A more detailed master chart of the computer program can 
be found in Appendix A. A detailed description of the optimization subprogram ‘ 
AESOP is presented in Appendix B. 
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OPTIMIZATION LOOP PERTURBS THE DESIGN VARIABLES 
TO OBTAIN THE GREATEST IMPROVEMENT IN RATING 


START 


FIXED ENGINE AND PROPELLANT DATA 
PLUS FREE DESIGN VARIABLES 


l‘ 


SPECIFIC COMBUSTOR 


CHARACTERISTICS 


CALCULATES PERFORMANCE, 
STABILITY, AND OTHER SPECI- 
FIC ENGINE CHARACTERISTICS 



USER-DEFINED WEIGHTING FACTOR 
CONSTANTS TO DEFINE AVERAGE 
COMBUSTOR CHARACTERISTICS 
REFLECTING THE DESIGNER'S 
VIEWPOINT ON THEIR REUTIVE 
SIGNIFICANCE 
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AVERAGE COMBUSTOR 
CHARACTERISTICS ‘ 

COMPUTES AVERAGE COMBUSTOR 
CHARACTERISTICS AS WEIGHTED 
AVERAGE OF THE APPROPRIATE 
SPECIFIC ENGINE CHARACTER. 



USER-DEFINED WEIGHTING FACTOR 
CONSTANTS FOR COMBUSTOR RATING 
BASED ON THE IMPACT OF AVERAGE 
COMBUSTOR CHARACTERISTICS ON 
VEHICLE PAYLOAD CAPABILITY/ AN 
IDEAL ENGINE WOULD HAVE A RATING 
OF ZERO 


COMBUSTOR RATING 

PROVIDES A SINGLE NUMERICAL 
MEASURE OF THE VEHICLE PAY- 
LOAD LOSS ASSOCIATED WITH A 
PARTICULAR DESIGN BY COMPU- 
TING THE SUM OF LOSSES 
ASSOCIATED WITH EACH WEIGH- 
TED AVERAGE COMBUSTOR 
CHARACTERISTIC 


[i 

OUTPUT OF OPTIMIZED COMBUSTOR 
DESIGN CHARACTERISTICS 


Figure l-'l. 


Conceptual Layout of the Optimization 
Computer Program 
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2 . INPUT DATA SETUP 


Figure 2-1 depicts the overall schematic of the AUTOCOM program with the 
respective input locations (heavy lines) . The input data setup is entered 
in nine blocks within the program (two data blocks are optional) . The data 
input sequence is presented in Figure 2-2, 

The input data setup describes a nominal engine design (for input data 
listing see Table 4-1, Section 4.1). There are the following nine input data 
types for the program: 

1. Informative data (title) 

2. Weighting factors 

3. Independent design variables 

4. Fixed (miscellaneous) combustor data 

5. Data for root guesses 

6. "Read** and **write" control 

7. Curve-fit data for propellant properties (optional) 

8. Optimization data 

9. Input for specific engine characteristics (optional) 

Provisions are made in the coding to consider 

Pump-fed or pressure-fed engines 

Engines with or without prebumers ; 

Various injector element types, e.g., coaxial or non-coaxial 

injection , showerhead 

Gaseous or liquid propellant state of propellants 

Hypergolic or non-hypergolic propellants 

Radial and/or circumferential baffles 

The weighting factors are normally set equal to "1**, and it is at the 
discretion of the designer to select his own inputs. (For the derivation of 
weighting factor constants see Section 5.5),. 
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INPUT CONSISTING OF 
FIVE CARD TITLES 


OPTIONAL 
DATA BLOCK 


OPTIONAL 
DATA BLOCK 



INPUT OF COMBUSTOR 
WEIGHTING FACTORS FOR 
SUBROUTINE CMWGTl 


INPUT OF DESIGN VARIABLE 
VALUES FOR SUBROUTINE 
CMDPVI 


INPUT OF MISCELLANEOUS 
COMBUSTOR DATA FOR SUB- 
ROUTINE CMMISC 


INPUT FOR ROOT SELECTING 
DATA FOR SUBROUTINE 
XROOTS 


INPUT FOR LOCAL CHARACTER- 
ISTIC READ AND WRITE CONTROL 
DATA FOR SUBROUTINE CMIRW 


INPUT FOR OPTIONAL CURVE- 
FIT DATA FOR PROPELLANT 
PROPERTIES FOR SUBROUTINE * 
FITCRV 


OPTIMIZATION DATA AND BOUNDARY 
VALUES FOR THE DESIGN VARIABLES 
FOR SUBROUTINE BAESOP 


OPTIONAL .INPUT FOR EACH CHARAC- 
TERISTIC COMPUTATION (SUB- 
ROUTINES SECFll THROUGH SECF82; 
ROOT MAPPING INSTRUCTIONS, ETC.), - 


Figure 2-2. Namelist Data Input Sequence 
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The program coding is arranged in such a way that the subroutines for computing 
specific engine characteristics (SECFll through SECF82) as well as for computing 
the combustor rating are only activated if the respective weighting factors 
(e.g., AFll, AF12, . . AFI , AFII, • . ) are not **zero.*' This permits a 

user controlled bypass of undesired specific engine characteristic computations 
by choosing the respective weighting factors to be "zero.** The average combustor 
characteristics modules for perfprm^ce, stability, etc. will be bypassed if 
no computation within their respective specific engine characteristics modules 
takes place. Figure 2-3 depicts the schematic of this feature and indicates the 
relevant weighting factor inputs. 

There are basically two options for applying AUTOCOM to an engine. These 
options, controlled by input data, are 

1. point design or analysis 

2. design optimization 

The difference between these options is in the number of iterations determined 
by data input. For details see Section 2.8. 

As mentioned before, it is possible to analyze selected engine characteristics 
by setting, for example, AF26 = 1, and bypassing all other engine characteristics 
(AFXX = 0) . 

There are provisions made in the AUTOCOM code for the input of characteristic 
data; for details see Section 2.9. 

The input data deck setup is presented in Figure 2-4. Details regarding the 
individual input steps are treated in the following sections. 


2.1 Title 

Corresponding to Figure 2-2, the **title** is the first place where data inputs 
have to be executed. A **title** consisting of five cards provides user selected 
information about 

1. The copputer run number and the engine type to be analyzed 

2. The type of run 

3. The thrust level 

4. The propellant combination 

5. Whether or not it is a point design or an optimization run 

There is no restriction on the type of information or on the number of card columns 
used. There must be five cards in the title; blanks are acceptable • A sanple is 
presented in Figure 2-5. 
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Specific Combustor 
Characteristics 



I 

I 



I 

I 



Figure 2-3. Logic of Bypass Arrangement for 
Combustor Characteristics and 
Ra- 






Figure 2-4. Deck Set Up with Namelist Names 



Figure 2-5. Sample Arrangement of Title Cards 
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2.2 Data Input Namelist COMWGT 


The equation to obtain average engine characteristics as well as the 
rating equation contain weighting factor constants whose values are selected 
by the desi^er. These weighting factors reflect his view regarding the 
importance of each characteristic. In the general program codings the . values 
are chosen to’,- be ”1” or depending on whether the subroutine is activated 
("1") or bypassed (”0") . In the following the specific engine characteristic 
weighting factor constants and their definitions are listed. (For derivation 
of the weighting factor constants for rating an engine ^ see Section 5.5). 


PROGRAM 

SYMBOL 

VALUE 

DESCRIPTION 

AFll 

Positive number 
or 0.0 

Weighting factor constant for fuel vaporization 
in the average performance characteristic 
equation. If a gas, then AF11=^0. 

AF12 

Positive number 
or 0,0 

Weighting factor constant for oxidizer vaporization 
in the average performance characteristic equation. 
If o, gaSj then AF12=0, 

AF13 

Positive number 
or 0.0 

Weighting factor constant for C*-efficiency in the 
average performance characteristic equation. 

AF20 

Positive number, 
or 0.0 

Weighting factor constant for the fuel side chugging 
in the average combustor stability characteristic 
equation . 

AF21 

Positive number 
or 0.0 

Weighting factor constant for the oxidizer side 
chugging in the average combustor stability charac- 
teristic equation. 

AF22 

Positive number 
or 0.0 

Weighting factor constant for pulsed operation in 
the average combustor stability characteristic 
equation. 

AF23 

Positive number 
or 0.0 

Weighting factor constant for non-pulsed operation 
in the average combustor stability charactersitic 
equation . 

AF24 

Positive number 
or 0.0 

Weighting factor constant for the fuel stability 
decay rate in the average combustor stability 
characteristic equation. 

AF25 

Positive number 
or 0.0 

Weighting factor constant for the oxidizer stability 
decay rate in the average combustor stability charac- 
teristic equation. 
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PROGRAM * 
SYMBOL 

. Value 

DESCRIPTION 

AF26 

Positive number 
dr 0.0 

Weighting factor constant for longitudinal inodes 
in the average combustor stability characteristic 

V ... 

^ . 

equation 

AF27 

Positive mimber 

Weighting factor constant for transverse modes 


or O.Q' - 

in the average combustor stability characteristic 



equation 

AF28 

Positive number 
or 0.0 

Weighting factor constant for combustion response 
and flow response stability in the average combustor 
stability characteristic equation 

AF29 

Positive number 
or 0.0 

Weighting factor constant for nonlinear, stability 
analysis in the average combustor stability charac- 
teristic equation. 

AF31 * 

Positive number 
or 0 . 0 " 

Weighting factor constant for fuel in the average 
pressure drop characteristic equation 

AF32 

Positive number 
or 0.0 

Weighting factor constant for oxidizer in the average 
pressure drop characteristic equation 

AF41 

Positive number 

Weighting factor constant for fuel and oxidizer 

- 

■ or 0 . 0 * ^ 

holes in the average injector con^ilexity character- 
istic equation. 

AF42 

Positive number 
or 0.0 

Weighting factor constant for the volume of the 
oxidizer dome in the average injector complexity 
characteristic equation. 

AF43 

Positive number 
or 0.0 

Weighting factor constant for the volume of the fuel 
dome in the average injector complexity character- 
istic equation. 

AF44 

Positive number 

Weighting factor constant for the length of the 


or 0. 0 

oxidizer holes in the average injector complexity 
characteristic equation. 

AF45 

Positive number 
or 0. 0 

Weighting factor constant for the length of the fuel 
holes in the average injector complexity characteristic 
equation . 

AF46 

0.0 

Weighting factor constant for the injector type in 
the average injector complexity characteristic 
equation. (At this time, the module SECF46 is not 


! 

activated. ) 
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1 

PROGRAM 

SYMBOL 

1 VALUE 

i DESCRIPTION 

AF51 

Positive number 
or 0.0 

Weighting factor constant in the average combustor 
length characteristic equation 

AF52 

Positive number 
or 0.0 

Weighting factor constant exponential in the average 
combustor length characteristic equation 

AF61 

• 

Positive number 
or 0 . 0 -N 

Weighting factor constant in the average combustor 
diameter characteristic equation 

AF71 

Positive number 
or 0.0 

Weighting factor constant in the average mixture 
ratio characteristic equation 

AF72 

i 

Positive number 
or 0.0 

Weighting factor constant exponehtial in the average 
mixture ratio characteristic equation 

AF81 

Positive number 
or 0. 0 

Weighting factor constant in the average baffle 
characteristic equation 

AF82 

Positive number 
or 0.0 

Weighting factor constant in the average baffle 
characteristic equation 

BF71 

> 

Greater than or 
Equal to 0/F 

Weighting factor constant: in the average 
characteristic equation 

AFI 

Positive number 
or 0.0 

Weighting factor constant of the performance charac- 
teristic in the rating equation 

AFII 

Positive number 
or 0. 0 

Weighting factor constant of the stability charac- 
teristic in the rating equation 

AFIII 

Positive number 
or 0.0 

Weighting factor constant of the pressure drop 
characteristic in the rating equation 

AFIV 1 

i 

Positive number 
or 0 . 0 

Weighting factor constant of the injector complexity 
characteristic in the rating equation 

AFV 1 

Positive number 1 
or 0 . 0 

Weighting factor constant of the chamber length 
characteristic in the rating equation 

AFVI 

Positive number 
or 0 . 0 

Weighting factor constant of the chamber diameter 
characteristic in the rating equation 

AFVII 

Positive number 
or 0.0 

Weighting factor constant of the mixture ratio 
characteristic in the rating equation 

AFVI I I 

Positive number 
or 0.0 

Weighting factor constant of the baffle characteristic 
in the rating equation 
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PROGRAM 

SYMBOL 

VALUE 

DESCRIPTION 

BFI 

Positive number 
or 0. 0 

Weighting factor constant exponential to the perfor- 
mance characteristic in the rating equation 

BFII 

Positive number 
or 0.0 

Weighting factor constant exponential to the stability 
characteristic in the rating equation 

BFIII 

Positive number 
or 0.0 

Weighting factor constant exponential to the pressure 
drop characteristic in the rating equation 

BFIV 

Positive number 
or 0. 0 

Weighting "factor constant exponential to the injector 
complexity characteristic in the rating equation 

BFV 

Positiye number 
or 0. 0 

i 

Weighting factor constant exponential to the chamber 
length characteristic in the rating equation 

BFVI 

Positive number 

Weighting factor constant exponential to the chamber 


or 0.0 

diameter characteristic in the rating equation 

BFVII 

Positive number 
or 0 . 0 

Weighting factor constant exponential to the mixture 
ratio characteristic in^the rating equation 

BFVI I I 

Positive number 
or 0.0 

Weighting factor constant exponential to the baffle 
characteristic in the rating equation 

CFII 

Positive number 
or 0.0 

Weighting factor constant exponential to the stability 
, characteristic" in the rating equation 


2.3 Data Input Namelist COMDPV 


Each specific engine characteristic is a function of the design variables 
entering into the combustor design. Basically, there are fourteen design 
variables. One of these is not activated at this time but is available for 
program extension. 

The combustor design variables are listed in the data input namelist 
COMDPV, consistent with the sequence coded in the program. Following this 
listing iS' the ^'Definition of circumferential and radial baffle arrangement,” 
Figure 2-6. 
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Figure 2-6. Definition of a Typical Circumferential 
and Radial Baffle Arrangement 



PROGRAM 

SYMBOL 

UNITS 

VALUE 

DESCRIPTION 

DF 

\ 

in. 

Positive 

number 

Minimum diameter of fuel orifices (deter- 
mines pressure drop of the fuel side of 
the injectoi); see Figure 2-7 ' * 

DX 

in. 

Positive 

number 

Minimum diameter of oxidizer orifices 
(determines pressure drop of the oxidizer 
side of the injector); see Figure 2-7 

NF 

- 

Positive 

integer 

Number of fuel orifices 

NOX 

- 

Positive 

integer 

Number of oxidizer orifices 

VOLF 

in^ 

Positive 

number 

Volume of fuel manifold 

VOLX 

in3 

Positive 

number 

Volume of oxidizer manifold • 

LF 

in. 

Positive 

number 

Length of fuel orifices; see Figure 2-7 

LOX 

in. 

Positive 
number 1 

Length of oxidizer orifices;, see Figure 2-7 _ 1 

- 

- 

--- 

Available for program extension 

LC 

* in . 

Positive 

number 

Length of combustion chamber (from injector 
face to throat area); see Figures 2-8 and 2-12 

DC 

in . 

Positive 

number 

Chamber diameter at injector face; see 
Figures 2-8 and 2-12 

OF 


Positive 

number 

Propellant mixture ratio (oxidizer/ fuel) 

BFLRAD 


Positive 
integer or 
0.0 

Number of radial baffles; .see Figure '2-6 

BFLCIR 


Positive 
integer or 
0.0- 

Number of circumferential baffles; see 
Figure 2-6 ' 
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2.4 'Data Input Namelist COMMIS 


The data ^input namelist COMMIS contains all miscellaneous inputs which 
define the combustor design as well as the propellant characteristics. The 
figures found subsequent to the COMMIS listing are as follows: 


Figure 

2-7 

Coaxial injector definition 

Figure 

2-8 

Combustor chamber geometry definition 

Figure 

2-9 

Matrix for the input selection of ETYP, ICOAX, 
and INJCMP 

Figure 

2-10 

. Input definitions for pump-fed engine with 
single preburner 

Figure 

2-11 

Input definitions for punQ)-fed engine with double 
preburners 

Figure 

2-12 

Preburner input parameter definitions 


PROGRAM 

SYMBOL 



DESCRIPTION 

ANGIMP 

degree 

0 to 180 

Impingement angle of injection 

CDF 

' ! v . ‘ 

Positive 

number 

Discharge coefficient for fuel orifices 
(product of the velocity coefficient and 
the contraction coefficient). For more 
about CDF see Section 4.1 

CDX 


Positive 

number 

Discharge coefficient for oxidizer orifices 
(product of velocity coefficient and the 
contraction coefficient). For more about 
CDX see Section 4.1 

CPF 

Btu/ 

lb-°F 

Positive 

number 

Specific heat^of fuel (condition at 
injector exit) 

CPX 

Btu/ 
lb-®F ' 

. Positive 
number 

Specific heat of oxidizer (condition at 
injector exit) 

CHIP 


Positive 

number 

Compressibility of the fuel, defined as 
the per cent change in volume produced 
by a per cent change in pressure, 

dv/V 

CHIP = dp/P 
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PROGRAM 

SYMBOL 


UNITS 


VALUE 


DESCRIPTION 


CHIX 

- 

Positive 

number 

Compressibility of the oxidizer, defined 
as the per cent change in volume produced 
by a per cent change in pressure 




PH TV - 

dp/P 

DFVEL 

in. 

Positive 

number 

Exit diameter of the fuel orifices (deter- 
mines fuel injection velocity]; see 
Figure 2-7 

DFIVEL 

in. 

Positive 

number 

Inner diameter of the fuel ring (serves 
to compute DFVEL); see Figure 2-7 

DF2VEL 

in. 

Positive 
number . 

Outer diameter of the fuel ring (serves 
to compute DFVEL) ; see Figure 2-7 

DTHRTI 

in. 

Positive 

number 

Combustion chamber throat diameter; see 
Figure 2-8 

DXVEL 

in. 

Positive 

number 

Exit diameter of the oxidizer orifices 
(determines oxidizer injection velocity); 
see Figure 2-7 

DXIVEL 

in. 

Positive 

number 

t. 

Inner diameter of the oxidizer ring (serves 
to compute- DXVEL) ; see Figure 2-7 

DX2VEL 

in. 

Positive 

number 

Outer diameter of the oxidizer ring 
(serves to compute DXVEL); see Figure 2-7 

D2 

in. 

Positive 

number 

Combustion chcimber diameter at the onset 
of the converging section; see Figures ' 
.‘2-8 and 2-12 

ETYP 


1,2,3, 
or 4 

Principal element type. (Also see Figure 2-9) 
1 = coaxial 
2. =/ unlike impinging 
. 3 =, like impinging 

4 = showerhead 

FLINE 

in . 

LF' 

Length of the fuel supply line. Ih the 
present computer program' the value 
for FLINE is the same as for LF 

FLAREA 

. 2 
in 

DF2.pi 

4 

Cross sectional area of fuel supply line. 
In the present computer program the value 
for FLXREA = (DF**2 *PI/4)NF 

F28TBS 

radians 

1.9 

Phase angle due to delay in burning 
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PROGRAM 

SYMBOL 

UNITS 

VALUE ; 

DESCRIPTION 

I COAX 

J 


0,1, or 2 

Coaxial injection indicator (also see 
Figure 2-9) ; 

0 = non- coaxial 

1 = outer ring is fuel 

2 = outer ring is oxidizer 

lEQF 

- 

1 or 2 

Fuel feed system indicator: 

1 = pump-fed system (constant flow) 

2 = pressure-fed system (constant pressure) 

lEQX 

- 

* 1 or 2 

Oxidizer feed system indicator: 

1 = pump- fed system (constant flow) 

2 = pressure-fed system (constant pressure) 

INJCMP 


' 0,1, 2, 3,4 
or 5 

Injector complexity index (see also Figure 
2-9) :. 

0 = showerhead injection with rigimesh 

1 = showerhead injection (fuel and oxidizer) 

2 = concentric tube injection (includes , 

coaxial impingement) 

3 = like impingement fuel (or oxidizer) 

with showerhead oxidizer (or fuel) 

4 = unlike triplet injection 
’ 5 = pentad injection 

ISTATF 

- 

1 or 2 

Fuel state at injector exit: 
r = gas 
2 = liquid 

ISTATX 

1 


1 or 2 

Oxidizer state at injector exit: 

1 = gas . ' ' j 

2 = liquid , 

ISTORF ! 


0 or 1 

Fuel hypergolic indicator: 

0 = non-hypergolic 

1 = hypergolic 

ISTORX 


0 or 1 

Oxidizer hypergolic indicator:' 

0 = non-hypergolic 

1 = hypergolic 

LAMBDF 

Btu/lb 

.Positive 
. number or 
0 

Latent heat of vaporization of fuel (for 
fuel condition at injector exit). If a » 

gas or above critical temperature, enter 
”0” 

LAMBDX 

Btu/lb 

Positive 
number or 
0 

Latent heat of vaporization of oxidizer 
(for oxidizer condition at injector exit) . 
If a gas or above critical temperature, 
enter "0” 








PROGRAM 

SYMBOL 

UNITS 

VALUE 

DESCRIPTION 

LG50F 

in. 

2.75 

1 

Generalized chamber length to vaporize 
50 per cent of fuel (see NASA TR R-67, 
Figure 28f) 

LG50X 

in. 

2.75 

Generalized chamber length to vaporize 

50 per cent of oxidizer (see NASA TR 1 

R-67, Figure 28f) 

LR 

1 


0 or 1 

Liner indicator in A. D. Little corre- 
lations (NASA CR- 72370) : 

1 = liner or absorber is present 
0 = liner or absorber is not present 

LI 

in. 

Positive 
number or. 
0 

Length of cylindrical section of the 
combustion chamber; see Figures 2-8 and 2-12 

MOLEF 

lb /mole 

Positive 

number 

Molecular weight of fuel (for definition 
of input values see Section 7.9.19) 

MOLEX 

Ib/mole 

Positive 

number 

Molecular weight of oxidizer (for defi- 
nition of input values, see Section 7.9.19) 

MUF 

Ib/in-sec 

Positive 
number or 
0 

Viscosity of liquid or gaseous fuel 
(condition at injector exit) . If a gas 
or above critical temperature, enter 

MUX 

Ib/in-sec 

Positive 
number or 
0 

Viscosity of liquid or gaseous oxidizer 
(condition at injector exit). If a gas 
or above critical temperature, enter ^*0*’ 

NTAUNl 

- 

Present 
value is 
.70 

Propellant interaction index. Used for 
modifying PINTER in subroutine TAUNl 

NTAUN3 

- 

Present 
value is 
.70 

Propellant interaction index. Used for 
modifying PINTER in subroutine TAUN3 

PCRITF 

Ib/in^ 

1 

Positive 
number or 
1.0 

Critical pressure of fuel. If a gas, 
PCRITF =1.0 

PCRITX 

i Ib/in^ 

Positive 
number or 
1.0 

Critical pressure of oxidizer. If a gas, 
PCRITX = 1.0 


20 







PROGRAM 

SYMBOL 


UNITS 


VALUE 


DESCRIPTION 


PTYP 


1,2, 3, 4, 
or 5 

* ' - i ’ ’ 

Propellant type indicator: 

1 = H2/02 

2 = RPl/LOX or JP5A/L0X 

3 = N2/H4 - UDMH/N204 

4 = Ethanol/LOX 

5 = Others 

RHOF 

Ib/in^ 

Positive 

number 

Density of fuel (condition at injectpr 
exit) 

NOTE: For an engine arrangement with 

prebumers , the density of the 
gaseous fuel/OX mixture at main 
chamber entrance becomes 




WXPRFR 

RHOX 

Ib/in^ 

Positive 

number 

Density of oxidizer (condition at injector 
exit) 

SIGLF 

Ib/in 

Positive 
number or 
0 

Surface tension of fuel (condition at 
injector exit). If a gas or above 
critical temperature, enter. "0" 

SIGLX 

lb /in 

! 

Positive 
number or 
0 

Surface tension of oxidizer (condition at 
injector exit). If a gas or above 
critical temperature, enter ’*0"' 

TCRF 

i 

°R 

Positive 
number or 
0 

Critical temperature of fuel. If a gas 
or above critical temperature, enter "0” 

TCRX ' 

°R 

Positive 
number or 
0 

Critical temperature of oxidizer. If a 
gas or above critical temperature, enter 
"0" 

THRUST 

Ibf 

• 

Positive 
number or 
0 

Effective thrust for respective engine 
operating conditions. If a prebumer is 
analyzed , enter "0" 

TLOF 

°R 

Positive 

number 

Droplet or gas initial temperature of 
fuel (condition at injector exit) . 

TLOX 


Positive 

number 

Droplet or gas initial temperature of 
of oxidizer (condition at injector exit) 

WFCOOL 

Ibm/sec 

Positive 
number or 
0 

Fuel mass flow rate for cooling (separ- 
ately injected into the combustion 
chamber upstream of the throat area) . 

If no cooling flow exists, enter "O’*. 

See Figures 2-10 and 2-11 for definition 
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Figure 2-8. Combustion Chamber Geometry Definitions 



Figure 2-9. Matrix for the Input Selection of ETYP, ICOAX, 
and INJCMP (crossed-out fields mean that no 
selection is possible) 

23 













wr ■ W7 ♦ wox 


r 



WT • WF ♦ WOX 






Figure 2-10. Input Definitions for 
Pump-Fed Engine with 
Single Preburner (Main 
thrust chamber and pre- 
burner are analyzed 
separately) 


Figure 2-11. Input Definitions for 
Piomp-Fed Engine with 
Double Preburners (Main 
thrust chamber and pre- 
burners are analyzed 
separately) 



Figure 2-12. Preburner Input Parameter Definitions 
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PROGRAM 

SYMBOL 

INPUT 

VALUE 

DESCRIPTION 

WFPREB 

Ibm/sec 

Positive 
number or 
0 

Fuel mass flow rate for prebumer. If 
no prebumer arrangement exists, enter 
"0”. See Figures 2-10 and 2-11 for definition 

WT- 

Ibm/sec 

Positive 

number 

Total propellant mass flow rate entering 
the combustion chamber. See Figures 2-10 
and 2-11 for definition 

WXPREB 

Ibm/sec 

) 

Positive 
number or 
0 

Oxidizer mass flow rate for prebumer. If 

no prebumer arrangement exists , enter 

”0". See Figures 2-10 and 2-11 for definition 

XLINE 

in. 

LOX 

Length of oxidizer supply line. In the 
present computer program the value for 
XLINE is the same as for LOX 

XLAREA 

in^ 

Nox 

4 

Cross sectional area of oxidizer Supply 
line. In the present computer program 
the value for XLAREA = (DX**2*PI/4)N0X 


2.5 Data Input Namelist NLROOT 

Subroutines SECF20, 21, 26, 27, and 28 compute low and high frequency 
combustion stability modes which require the solution of complex characteristic , 
equation roots. The root finding procedure is assisted by either nominal root 
guesses stored in the program code or by root guesses specified by the designer as 
inputs in input namelist NLROOT. Root guesses for namelist NLROOT are convenient 
for repetitive runs. Previously found roots are used to facilitate the 
searching process; i.e., to save computer time. In addition, various kinds 
of control statements are inputs in NLROOT. 

. There are three types of control statements; namely 

1. general control statements . 

2. control statements for individual subroutines - 

3. control statements for root guesses 

These control statements are discussed in the sections to follow. 











2.5.1 General Control Statements. The general statements 


ROOTS = .FALSE. 

ROOTS = .TRUE. • 

indicate if there. are any root guesses stored on tape. The general statements 

IPROOT = . FALSE . 

IPROOT = .TRUE. 

are set true if a printout of the root table summary is desired. 

2.5.2 Control Statements for Individual Subroutines . To determine if root 
guesses for the individual subroutines are stored in the program code (.FALSE.) 
or if root guesses are used in input data (.TRUE.), the respective, control 
statements become 

- IGRF20 = .FALSE., (or .TRUE.) 

IGRF21 = .FALSE, (or .TRUE.) 

IGRF26 = .FALSE, (or .TRUE.) 

IGRF27 = .FALSE, (or .TRUE.) 

IGRF28 = .FALSE, (or .TRUE.) 

The maximum number of roots possible to be solved is indicated as follows: 

NRF20 =20 

NRF21 = 20 . , . 

NRF26 =20 

NRF27 = 10 , 

NRF28 =10 

It is at the user's discretion to select any smaller number of roots. 

The maximum number of trials possible for each root is indicated as 
follows: 

NTF20 = NRFXX * NTFXX <: 20 
NTF21 = NRFXX * NTFXX < 20 
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/ 


NTF26'= NRFXX * NTFXX < 20 . 

NTF27 = NRFXX * NTFXX <20 
NTF28 = NRFXX * NTFXX < 20 ' 

Here,' too, it is up to the judgment of the user, possibly to reduce the number 
of trials for a given root. 

The search range for decay rate in the root finding is indicated by the 
statements 

' ALRN27 = 30 ■ - * ■ ^ 

" . ' ' ALRN28 = 30: ’ . ; . . . 

Other range definitions (foi SECF20, 21, "26) are built in the program code and; 
therefore, cannot be varied. 

In subroutine SECF27 the maximum number of modes of the oscillation 
constant, is coded in the program; therefore, 

NF27 =10 

In subroutine SECF28 the maximum number of Bessel argiments, m, is coded in the 
program; therefore. 


NB28 = 10 ■ ‘ 

Note that if baffles are' present the program will automatically determine NF27 
and NB28. 


2.5.3 Control Statements for Root Guesses . A typical root guess in data 
input NLROOT has the form 


G R 

\ / 

V 

t 

guess root 


F 27 



subroutine SECF27 


^ R 


( 1 , 2 ) 


real part 
of the root 


first 

frequency 


second 

root 
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Corresponding 


GRF27I(1, 2) = 


i 

imaginary part 
of the root 



The following Tables 2-1 to 2-4 list the root guesses for subroutines 
SECF20, 21, 26, 27, and 28, respectively, stored in the program code as, well 
as the sequence for the data input namelist NLROOT. .. 


2.6 Data Input Namelist COMIRW 

During the construction of the computer program it became necessary to* 
have the option ^writing out" (IRWXX = 1) the results of a selected subroutine 
(DBFXX or BPAL namelist), to "read in" (IRWXX = 2) specific data; or by setting 
IRWXX = *0 to ignore "writing out" and "read in". Subroutines SECFll through 
SECF82 can be controlled in this way as well as several other subroutines. 

Table 2-5 depicts the options for IRWXX contained in the program code. For 
more details about "writing out" and "read in" see Section 2.9. 

2.7 Data Input Namelist CRVFIT 

The curve fitting capability in the computer code is necessary because 
the four: parameters 

Isp = specific impulse 

K = specific heat ratio 

M = molecular weight of exhaust gas 

T^ = temperature in combustion chamber 

are, for optimization purposes, quantities as a function of the mixture ratio 
0/F. For point design, however, the listed parameters are fixed input values. 
Therefore, two possibilities for data input are considered: 

1. Use of data stored in program. The input is 

$CRVFIT 

$ 

2. Use of fixed input values for point design if 
the relationship of the four propellant param- . 
eters as a function of' 0/F is not available. 

The input cards are then 
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Subroutine Subroutine I Subroutine Subroutine 

SECF26 (Oxidizer) - SECF26 (Fuel) ■ SECF21 (Oxidizer) SECF20 (Fuel) 



REAL 

imaginary 

^AME 

IN 

PROGRAM 

root guess 
STORED IN 
PROGRAM 

NAME 

IN 

PROGRAM 

ROOT GUESS 
STORED IN 
PROGRAM* 


GRF20R(1) 

GRF20R(2) 

GRF20R(3) 

GRF20R(4) 

GRF20R(5) 


GRF21R(11) 

GRF21R(12) 

GRF21R(13) 

GRF21R(14) 

GRF21R(15) 



GRF20I(1) 

0.0 

GRF20I (2) 

-tt/t 

GRF20I(3) 

7T/t 

GRF20I(4) 

to 

GRF20IC5) 

-0) 

GRF21I(11) 

0.0 

GRF21I(12) 

-tt/t' 

GRF21I(13) 

tt/t* 

GRF21I(14) 

(0 

GRF21I(15) 

-0) 


GRF26R(1) 

-O.OIOI 

GRF26I(1) 

PIFAC 

GRF26R(2) 

-0.0101 

GRF26I(2) 

-PIFAC 

GRF26R(3) 

-0.0101 

GRF26I(3) 

PIFAC 

GRF26R(4) 

-0.0101 

GRF26I(4) 

-PIFAC 

GRF26R(11) 

0.0101 

GRF26I(11) 

PIFAC 

GRF26R(12) 

0,0101 

GRF26I(12) 

-PIFAC 

GRF26R(13) 

0.0101 

GRF26I(13) 

PIFAC 

GRF26R(14) 

0.0101 

GRF26I(14) 

-PIFAC ‘ 
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* K|. ci random nui:^er between 0 and *0.001 
*♦ ZUir « u)j/2Tr 


TABLE 2-3. Root Guesses for Subroutine SECF28 



























































SNLROOT 

roots=.true:,* 

ROOtS=. FALSE.* 

IPROOT=.TRUE.» 

NRKPOslfNTFaOsa.NRFZlsatNTFaisa* 

NRF26=5*NTF26=5. 

NRF?6=3»NTF26=3* 

ICRF28=2.NTF2B=2. 

NF27=7* Ng?8=5* 

NRF27=3.NTF27=3f 
IOPF20 =. false.* 

IGRF20S.TPUE.* 

IGPF21=.FALSE.» 

ir,PF21 = .TRU£.*' . - , 

IGRF26=. False.# r . - 

IGRF26=.TWUE.. 

IGRF27=.TRL’E.* 

1GRF27=. FALSE.* 

IGPF2«=. FALSE.* 

IGRF28=.TRUE.* 

ALRN28=lO.O* ALPN27=15.0* 

ALRN28=10.0* ALRri27=100.0* 

GRF20R(1) = .16369# .16399# 

GRF?Oin) = -.A3A81, .43/^?.!* 

GRF21R(1) = 2.026219E-02* 2.026219E-02# -.850782* 
GRF21K1) = ,833296# -.833296* 2.48576# 

GRF26R(1) = -.37924# 

6RF26IU) = -7.83507E-08* 

GRF26R(21 = -.33717* 

GRF26K2) = 3.2757# 

6RF26Rt3> = -.37969# 

GRF26K3) = 1.27032E-07* 

GRF26R(4) = -.33456# 

GRF26IC4) = 3.2736* 

GRF26R(11) = -.37924# 

GRF26I(ll) = -7.83507E-08* 

GRF26R(12) = -.33717# 

GRF26Itl2) = 3.2757# 

6RF26Ra3) = -.37969* 

GRF26K13) = 1.27032E-07* 

GRF26R(14) = -.33456* 

GRF261il4) s 3.2766* 

GRF28R(l*l) = 1.61195# 

6RF281 (1 .1) = -.145074# 

GRF2«R(1#2) s -1.61195* 

GRF28I (1*2) = -.145074* 

GRF28R(2*1) = 2.66034, 

GRF28K2*!) = -.229889* 

CRF28R(2#2) = -2.66034* 

GRF28I(2*2) = -.229889* 

GRF28R(3*1) = 3.66597* 

GPF28I (3.1) = -.43242# 

GRF28R(3#2) = -3.66597# 

GPF23I (3#2) = -.43242# 

GRF28R(4*1) = 3.33687# 

GRF28I (4.1) = -.34992# 

GPF28R(4#2) = -3.33607# 

GRF2RI (4*2) = -.34992# 

GPF28R(5#1) = 6.783433E-02* 

G«F28I(5^1) = -6.872498E-02# 

GRF28RC5.2) = -6.78J433E-02# 

GPF28K5.2) = -6.872498E-02# 


TABLE 2-4. Sample for Control Cards for Data 
Input Namelist NLROOT 
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IRWXX SUBROUTINE NAMELIST 


OPTIONS 

IRWACW = 0, I 
IRWFLO = 0, I 
IRWFRP = 0, I, 2 
IRWL50 = 0, I 
IRWNLS =0, I, 2 
IRWPAL = 0, 1 
IRWTAU = 0, I 
IRWVAP =0, I, 2 
IRWII =0, 1, 2 

IRW82 = 0, I, 2 


NAME NAME 

ACWAVE DBACW 

FLORES DBFLO 

FRPROP DBFRP 

LFIFTY DBL50 

NLSTAB DBNLS 

POSTAL BPAL 

TAUN DBTAU 

VAPRIZ DBVAP 

SECFIl DBFIl 



SECF82 I DBF82 


TABLE 2^5. Options for IRWXX Control Switch Contained in the 
Program Code 
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$CRVFIT 

NPISP=0 

NPK=0, 

NPMOLW=0, ' ' 

NPT1=0 
PCISP(1)=. . 

PCK(l) = . . 

PCM0LW(1)=, . 

PCT1(1)=. . 

$ 

The program symbols ‘are described in the list at the end of ■Secti6h'’2.'7;. 

2.7.1 How to Store New Propellant Data . If optimization runs have to be 
perfoiined, the relationship of the four quantities of Isp/ sxid as' 

function of the mixture ratio 0/F has to be stored in the program. For this * 
purpose, subprogram FITCRV uses four pertinent propellant polynomial 
expressions. The independent variable of these polynomial functions is the i 
mixture ratio 0/F. (For more details see Section 6.6.3). The polynomial 
coefficients must be supplied by the designer either directly or, indirectly 
by specification of the curves' in tabular form. For a new propellant combi-- 
nation, the sample data deck is as follows: 

$CRVFIT 

XFIT=0F2, OF2 . . . OFjyjpTg 

YFIT=ISPi, ISP2 . . . ISPnpts 
NPTS=10, 

KP0LY=2, 

KIND=1, 

ICRVFT=1, 

LOOP=.TRUE. ^ ^ 

$ ■ ' 

$CRVFIT • ‘ ? ■ ' ' ' 

XFIT=OFi, OF2 . . . OFnpTS 
YFIT=Ki, K2 . . . KnPTS 

NPTS=10, ' ^ 

KP0LY=2, 

KIND=1, 

ICRVFT=2, 

LOOP=.TRUE. 

$ ^ ' . 

$CRVFIT 

XFIT=OFi, OF2 . . . OFnptS 
YFIT=Mi, M2 . . . MnptS 
NPTS=10, 

KP0LY=2,- 

KIND=1, > : 

ICRVFT=3, 

LOOP = .TRUE. ‘‘ ' 

; Continued 
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$CRVFIT 

XFIT=OF,,OF„ OFnptq 

YFIT-n ,Tl| TiKI 

NPTS=10, 

KPOLY=2, 

KIND=1, 

ICRVFT=4, 

LOOP= . FALSE . 

$ 


program" 

SYMBOL . 

. DESCRIPTION 

ICRVFT 

=1, fit I._ to 0/F; =2, fit k to 0/F; 

=3, fit M ^ to 0/F; . ,=4, fit TI to 0/F . 

lOUT 

Output; polynomial values at input points . » 

ipnaml/, 

) 

Causes to print out the namelist, . 

JPOLY 

Number of leading polynomial zero coefficients 

KIND 

\ 

Kind of fit : 

= 1 for LINX, LINY 
= 2 for LOXG, LINY 
= 3 for LINX, LOGY 
= 4 for LOGX, LOGY 

KPOLY 

Degree of approximating polynomial 

LOOP 

Logical variable indicating to loop back for the 
next case 

LWGHT 

Weight switch : 

= 0 for WEIGHT (I) =1.0 
= 1 for WEIGHT (I) = DATA 

NPISP 

Polynomial degree of the specific impulse, , as a 

function of the mixture ratio 0/F ^ 

NPK 

Polynomial degree of the, specific heat ratio, k*, as a 
function of the mixture ratio 0/F 

NPMOLW 

Polynomial degree of the molecular weight of the combustion 
gas, M, as a function of the mixture ratio 0/F. 

NPTS 

Number of points on original curve 

NPTl 

Polynomial degree of the combustion temperature, Tl, 
as a function of the mixture ratio 0/F 


(continued) 
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PROGRAM 

SYMBOL 

DESCRIPTION 

PCISP 

Input polynomial coefficient for specific 
impulse, I^^ 

PCK 

Input polynomial coefficient for specific 
heat ratio, k 

PCMOLW 

Input polynomial coefficient for the molecular; 
weight of the combustion gas, M 

PCTl 

Input polynomial coefficient for combustion 
temperature, T1 

WEIGHT 

Weighting matrix to be applied to the raw data 

XFIT 

X values 

YFIT 

Y values • ' . 'i 


2.8 Data Input Namelist lAESOP 


In the present program code there are 14 design variables which are input 
data in namelist COMDPV. They are set in subroutine CMBDPV in the following 
sequence: . . 


DF- 
DX 
NF . 

NOX 
VOLF 
VO LX 
LF 
LOX 

(available) 

LC 

DC 

OF 

BFLRAD 

BFLCIR 


ALPHA (1) 

" ( 2 ) . 

" (- 3 ) 

" (4) 

. " (5) 

" ( 6 ) 

" (7) 

" ( 8 ) 

" (9) 

” (10) 

" ( 11 ) 

V. (12) 

" (13) 

" (14) 


For a deck setup, the boundaries (low and high) for these variables have to 
be determined and are inputs in lAESOP in the following manner: 


ALPHA(l) = ^ (14) 

ALPLO(I) = ^ (14) 

ALPHI(l) = ^ (14) 
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In addition, the number of iterations for the optimization loop has to 
be an input such as 


MAXJJJ = 10 

for ten iterations or in case of a point design 

MAXJJJ = 0 

Table 2-6 depicts an input sample for lAESOP. A description of the program 
symbols is given in the following list. More details about the optimization 
program AESOP and its . application are presented in Appendix B. 


$ I AESOP 

nalpha=ia, 

NFUNC = 15» 

IRFPET=20» 

LUFIN=2» ITFIN=0» 

NUM0PT=2» METHOP(1)=10»2. 

MAXRRS=5» 

■ MAXRRS=A» 

IPGAIN=1» 

IPSRCH=1» IPACYC=1, IPFCYC=1» 

MAXJJJ = 1» 

MAXJJJ = 2« 

MAXJJJ = 0» 

ALPL0(1)= .020* .020» 100. » 100. » 100000. » 100000. » .050t .050* 
ALPL0(9)= 0.* 10. ♦ 70.» 2.398» ?.♦ l.» 

ALPHA(1)= .195» .242. 2000.* 2000.* 100000.* 100000.* 2.0* 5.0* 
ALPHA(9)= 0.» 50.* 94.* 2.398* 7.* l.» 

ALPHI(1)= .50* 1*0* lO.OE+06* lO.OE+06* 100000.* 100000.* 12.* 18. 
ALPHI(9)=0.* 150.* 140.* 2.398* 7.* 1.* 

$ 


TABLE 2-6. Input Sample for lAESOP 


PROGRAM SYMBOL 

DESCRIPTION 

ALPHAi 

The nominal parameter values 

ALPHI^ 

Upper bounds on each parameter search range 

ALPLOi 

Lower bounds on each parameter search range . 

IPACYC 

Will cause print out of the control parameter values 
at the end of each optimization cycle 

IPFCYC 

Will cause print out of the function values at the 

end of each optimization cycle . .. 

^ ^ (continued) 
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PROGRAM SYMBOL 


DESCRIPTION 


IPGAIN 

Will cause print out of the function and control parameter 
values after each evaluation in which there was a gain 
in performance 

IPSRCH 

N 

Will cause print out of the function and control param- 
eter values at the end of each search 

IREPET 

The maximum number of times the search sequence (optimi- 
zation cycle) defined in METHOPj^ will be utilized^ 

ITFIN 

Controls content of final summary report 

LUPIN 

Tape or disk unit used to store the convergence histories 
of the function and control parameter values for tabular 
listing at the end of the optimization process . 

MAXJJJ 

The maximum number of iterations 

^ maxr'rs 

The maximum number of random rays, one or, two sided, 
investigated each time the optimization cycle 
requests a random ray search 

METHDP . 
1 

The search sequence by numeric identification 

NALPHA 

The number of parameters available for optimization 

NFUNC 

The total number of functions being defined and 
computed in the system model - ^ 

NUMOPT 

The number of optimization techniques to be employed 


2,9 Data Input Namelist for Characteristic Data 

The use of the data input namelist for characteristic data allows the 
designer to treat several program modules essentially in an independent 
fashion. By this it is possible to investigate the influence of selected 
parameters on the results of the respective module. Section 2.6 described 
the method for using the control switch for "reading in" characteristic data 
by setting IRWXX = 2. By doing this the data can be read in as follows: 

$ ‘ ■ ' . ■ : ^ 

DBXX 


s A list of each DBXX namelist available in the program with a description; 
of the namelist variables follows: ~ ^ . 
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CONTROL 

NAMELIST SWITCH FOR SUBROUTINE NAMELIST 

NAME READ AND NAME VARIABLES 
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3. PROGRAM OUTPUT 


The present program coding provides the print out of seven summary 
tables with parameter descriptions controlled by the following subroutines 


1. 

Subroutine 
Table 3-1 

OUTAEC, 

2. 

Subroutine 
Table 3-2 

OUTFDV, 

3. 

Subroutine 

3-3 

OUTGEP, 

4. 

Subroutine 
Table 3-4 

OUTSCC 

5. 

Subroutine 

OUTSRS 

6. 

Subroutine OUTWAC 
Characteristics ;** 

7. 

Subroutine 
Table 3-7 

OUTWTR 


"Average Engine Characteristics 
"Final Values of the Design Variables,"' 
"General Engine Parameters," Table 
"Specific Combustor Characteristics," 

"Stability Root Summary,’-* Table 3-5 

**Weighting Factors in Average 
Table 3-6 - - 

**Weighting ^Factors in Rating Equation,** , 


3.1 Data to Regulate Program Output 


The print out of the seven summary tables, is automatically .regulated 
with each run. Other informative print outs are activated when using the 
control switch for **read in** and **print out** IRWXX = 1 or 2 » (see Sections 
2.6 and 2.9). 



KUMRER ^yy 


engine type no,! 

TEST CASE 

thrust s 15000 POUND 
0? / H2 PROPELLANT 
CHECK OUT RUN 

AVERAGE ENGINE CHARACTERISTICS 

performance characteristic 
stability characteristic 

PRESSURE DROP CHARACTERISTIC 
INJECTOR COMPLEXITY CHARACTERISTIC 

length characteristic 
chamber diameter characteristic 
chamber mixture ratio characteristic 
combustor baffle characteristic 


»«««»««««««««««««««« 

TABLE 3-1. Typical Printout for Average Engine 
Characteristic^' (Subroutine OUTAEC) 


NUMBER AYY ENGINE TYPE NO.l 

TEST CASE 

THRUST = ISOOO POUND 
02 / H2 PROPELLANT 
CHECK OUT RUN 

final values of THE DESIGN VARIABLES 


DIAMETER OF FUEL ORIFICES 

OF 

ae 

,129200 

DIAMETER OF OXIDIZER ORFICES 

ox 

s 

fl,AOOOOOE-02 

NUMBER OF FUEL ORIFICES 

nf 

s 

216,000 

NUMBER OF OXIDIZER ORFICES 

nox 

8 

216.000 

VOLUME OF FUEL MANIFOLD 

VOLF 

3 

10.0000 

VOLUME OF OXIDIZER MANIFOLD 

VOLX 

3 

10.0000 • 

length of fuel ORIFICES 

Lf 

8 

6.000000E-02 

length OF OXIDIZER ORFICES 

.LOX 

8 

.400000 

element TYPE 

etyp 

X 

1.00000 

length of chamber 

LC 

s 

11.2200 

chamber diameter 

DC 

3 

10.2800 

mixture RATIO (OXIDIZER/FUEL) 

OF 

8 

5.05990 

RADIAL baffles 

bflrao 

8 

0. 

circumferential BAFELES : 

BFLCIR 

3 

0. 


FINAL COMBUSTOR RATING « 391.920 


TABLE 3-2, Typical Printout for Final Values of 

the Design Variables (Subroutine OUTFDV) 


rating component 


FI 

= 

.463970 

69,5955 

FII 


-37.4853 

4.71049 

Fill 

= 

71,4208 

1.49876 

FIV 

= 

8.87860 

3.11377 

FV 

r 

2.18268 

173.002 

FVI 

= 

2.00000 

lAO.OOO 

FVII 

= 

4.000000E-03 

1.729552E’ 

FVIII 

s 

0. 

0. 


«««-»« 


COMBUSTOR RATING = 391.920 
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NUMBER 4YY 


ENGINE type NO.l 
TEST CASE 

THRUST = 15000 POUND 
02 / H2 PROPELLANT 
CHECK OUT, RUN 

GENERAL ENGINE PARAMETERS 


COMPUSTOR thrust FORCE 


THRUST 

= 

15000.0 

propellant specific impulse 


ISP 

= 

AAA. 840 ‘ 

chamber pressure at injector head 


PCI 


372.817 

TOTAL propellant FLOW RATE 


WT 

= 

33.7200 

FUEL weight Flow rate’ 


WF 


5.56AA5 

OXIDIZER weight FLOW RATE 


WOX 


28.1556 

combustion temperature in chamber 


TCOMB 

= 

S76A.91 

RATIO OF SPECIFIC HEAT OF COMBUSTION GAS 


SPHEAT 

= 

1 .24900 

GAS constant OF COMBUSTION GAS 


R6AS 


127.131 

IDEAL thrust COEFFICIENT . 


CFIDEL 

= 

1.939CI 

acoustical length OF chamber 


LOUM 


8.81574 

molecular weight of COMBUSTION GAS 


MOLWT 


12*1450 

MEAN RESIDENCE TIME OF GAS IN CHAMBER 


THETAG 

= 

8.568235E 

SPEED OF SOUND IN CHAMBER 


CS 

= 

65123.3 

COMBUSTION chamber MACH NUMBER 


MC 

= 

.201078 

INJECTOR PRESSURE DROP FOR FUEL 


DELPF 


82.9049 

INJECTOR PRESSURE DROP FOR OXIDIZER 


DELPX 

= 

48.4525 

COMBUSTION chamber VOLUME 


VC 


681.565 

average VELOCITY OF GASES IN CHAMBER 


VELC 


13094.9 

FUEL INJECTION VELOCITY 


VFUEL 

= 

13013.0 

OXIDIZER INJECTION VELOCITY 


VOXIO 


339.511 

TOTAL AREA OF FUEL INJECTOR ORIFICES 


TAF 

= 

2*83184 

TOTAL AREA OF OXIDIZER INJECTOR ORIFICES 


'TAX 

= ■ 

1.19702 

chamber length TO VAPORIZE 50 PER-CENT OF 

FUEL 

L5FUEL 


0. 

chamber LENGTH TO VAPORIZE 50 PER-CENT- OF 

OXIDIZER 

L50XID 

= 

.649491 

baffle LENGTH BASED ON LARGEST L50 


8AFLEN 

3 

1.29898 


TABLE 3«3. Typical Printout for General Engine Parameters 
(Subroutine OUTGEP) 


N'UMf'ER 4YY 


K'^GINE • TYPE NO.l 
TEST C&SE 

THPuST = 15000 POUND 
02 / H? PROPELLANT 
CHECK OUT RUN 


SPECIFIC COMBUSTOR CHARACTERISTICS 




PER CENT MASS FUEL VAPORIZED 

FU 

=. 

0. 

PER CENT mass of OXIDIZER . VAPORIZED 

F12 

= 

99.8845 

C* EFFICIENCY MIXING MODEL 

F13 

= 

100.000 

efficiency pulsed COMBUSTORS 

F14 

= 

91,3821 

c* efficiency non-pulsed combustors 

F15 

= 

73.7667 

fuel system chugging decay rate 

' F20 

= 

-13004.4 

OZIOIZER system chugging DECAY PATE 

F21 

= 

-259.892 

PULSED instability CHARACTERISTIC 

F22 

= 

-37.4853 

non-pulseo instability characteristic 

F23 

= 

0. 

dykema fuel stability decay rate 

F24 

= 

0. 

dykeha oxidizer stability decay rate 

F25 

= 

-99.9432 

STABILITY longitudinal TIME LAG 

F26 

= 

-6182.25 

stability transverse TIME LAG •. 

F27 


-413.671 

STABILITY LRC RESPONSE FUNCTION 

F28 

= 

-356.492 

STABILITY PRIEM LINEAR ANALYSIS 

F29 

= 

-160.494 

FUEL PRESSURE DROP CHARACTERISTIC 

F31 


82.9049 

OXIDIZER PRESSURE DROP CHARACTERISTIC 

F32 

= 

48.4525 

FtiEL PLUS OXIDIZER HOLES CHARACTERISTIC 

F41 


431.000 

OXIDIZER DOME VOLUME CHARACTERISTIC 

F42 


10.0000 

fuel DOME VOLUME CHARACTERISTIC 

F43 . . 

■ = 

10.0000 

OXIDIZER HOLE LENGTH CHARACTERISTIC 

F44 ' 

= 

.400000 

fuel HOLE LENGTH CHARACTERISTIC 

F45 

■ = 

6.0000^00E 

INJECTOR TYPE COMPLEXITY CHARACTERISTIC 

. . F46 


0.- 

INJECTOR LENGTH CHARACTERISTIC 

F51 

= 

2.18268 

CHAMBER DIAMETER CHARACTERISTIC 

F61 

' = 

2.00000 

mixture RATIO CHARACTERISTIC 

F71 

= 

5.05990 . 

RADIAL baffles 

F81 

. = 

0. 

circumferential BAFFLES 

F82 

= 

0. 


TABLE 3-4, Typical Printout for Specific Combustor Characteristics 
(Subroutine OUTSCC) 
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TABLE 3-5. Typical Printout for Stability Root Summary 
(Subroutine OUTSRS) 


NUMBER 4YY 


engine type NO.l 
TEST Case 

THRUST = ISOOO POUND 
02 / H2 PROPELLANT 
CHECK OUT RUN 


WEIGHTING TACTORS IN AVERAGE CHARACTERISTICS 


CONSTANTS 

IN 

PERFORMANCE CHARACTERISTIC 

AFll 

s 

0. 




- AF12 

S 

l.OOOOO 




AF13 

= 

1.00000 




AF1<* 

= 

l .OOOOOOE- 




AF15 

= 

I .OOOOOOE< 

CONSTANTS 

IN 

ENGINE STABILITY CHARACTERISTICS 

AF20 

= 

1.00000 




AF21 

's 

1.00000 




AF22 

= 

uooooo 




AF23 

= 

0. 




AF2A 

_ =‘ 

0. 




AF25 

' = 

1.00000 




AF26 

_= 

1.00000 




AF27 

S. 

1.00000 




AF2A 

= 

1.00000 




AF2R 

= 

1 .00000 

CONSTANTS 

IN 

PRESSURE DROP CHARACTERISTICS 

AF31 

= 

2.00000 




AF32 

s: 

1.00000 

CONSTANTS 

IN 

THE INJECTOR COMPLEXITY CHARACTERISTICS 

AF41 

5 

R.OOOOOOE' 




AF4? 

= 

.125000 




AF43 

= 

,125000 




AF44 

= 

5.62500 

CONSTANTS 

IN 

THE ENGINE LENGTH CHARACTERISTICS 

’ AF51 

= 

1.00000 




AF5? 

= 

1 .00000 

CONSTANTS 

IN 

THE COMBUSTOR DIAMETER CHARACTERISTICS 

AF6I 

S 

1.00000 




AF6? 

= 

1.00000 




AF63 

= 

0. ■ ^ 

CONSTANTS 

IN 

THE MIXTURE RATIO CHARACTERISTICS 

AF71 

s 

4.00000 




AF72 

= 

2.00000 

CONSTAfiTS 

IN 

THE BAFFEL CHARACTERISTICS 

AF81 

ts 

0. . ' 




AF82 

= 

0. 

TABLE 

3- 

6. Typical Printout for Weighting 

Factors 

in 



Average Characteristics (Subroutine OUTWAC) 


NUMBER 4YY ENGINE TYPE NO.l 

TEST CASE 

THRUST = 15000 POUND 
0? / H2 PROPELLANT 
CHECK OUT RUN 

WEIGHTING factors IN RATING EQUATION 


‘V 


EXPONENTIAL ON PERFORMANCE CHARACTERISTICS 

AFI 

= 

15C.OOO 

CONSTANTS 

OF 

STABILITY CHARACTERISTICS 

AFII 


200.000 






BFII 

s 

0. 






Cfll 

= 

.100000 

CONSTANTS 

OF 

THE 

PRESSURE 

DROP CHARACTERISTIC 

AFIlI 


6.900000E' 






BFIII 

= 

1 .80000 

CONSTANTS 

OF 

THE 

INJECTOR 

complexity characteristic 

AFIV 

=: 

3.<950000E' 






BFIV 

= 

2,00000 

CONSTANTS 

OF 

the 

CHAMBER 

LENGTH characteristic 

AFV 


40.5000 






BFV. 


1.86000 

CONSTANTS 

OF 

the 

CHAMBER 

diameter characteristic 

AFV I 

= 

. 35.0000 






BFVI 

= 

2.00000 

constants 

OF 

THE 

MIXTURE 

RATIO characteristic 

AFVr I 


66.0000- 






BFVII 

3 

1 .16000 

CONSTANTS 

IN 

THE 

BAFFEL CHARACTERISTICS 

AFVIII 

3 

0. 






BFVIII 

~ 

0. 

TABLE 

3- 

7. 

Typical 

Printout for Weighting 

Factors 

in 



Rating Equation (Subroutine OUTWTR) 
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4. ENGINE. TEST CASE 


NOTE : 

I The following engine test run has been performed 
in the later part of 1971. Program modifications 
have been implemented since. Therefore, the print- 
outs do not show these modifications; however, the 
numerical results would still be unchanged with the 
exception of Subroutines SECF14 and SECF15 (computes 
c*-efficiency based on A. D. Little, Reference 7-4) 
which have been eliminated in the new AUTOCOM 
version at the request of NASA's technical monitor. 


4.1 Data Input for Test Case 


A 15 ,000-lbf-thrust LH 2 /LO 2 pump-fed rocket engine has been selected 
for demonstration of the AUTOCOM program. Combustor design variables for 
this engine (with no baffles present) are: 


1. Fuel Orifice' Diameter 


2. OX Orifice AP Diameter 

3. Number of Fuel Elements 

4. Number of OX Elements 

5. Vol. of Fuel Manifold 

6. Vol. of OX Manifold 

7. Length of Fuel Orifices 

8. Length of OX Orifices 

9. Length of Combustion Chamber 

10. Diameter of Combustion Chamber 

11. Propellant Mixture Ratio 


.195 inches O.D., .145 inches 

I. D. (.129 inches diameter 
equivalent hole) 

.084 inches 
216 
216 

10 in3 
10 in3 
.06 inches 
.40 inches 

II. 22 inches ■ 

10.28 inches 

5.06 Ibm OX/lbm Fuel 


Other pertinent but fixed design parameters include: 


Propellants 
Element Type 

OX Orifice Velocity Diameter 
Element Impingement Angle 
Total Propellant Flow 
Thrust 

LOX Temperature 


Hydrogen and LOX 

Concentric Tube (Hydrogen on outside) 

.11 inches 

0 degrees 

33.72 Ibm/sec 

15000 Ibf 

Boiling 
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H 2 Temperature 349® R 

Throat Diameter 5.14 inches 

This engine runs at the following measured conditions: 


83. 1 psi 

48. 3 psi f 

98,6 per cent ' . 

444 Ibf-sec/lbm . . 

396.4 psi ' 

These running conditions were used to check the nominal engine description 
in the combustor synthesis and to determine the discharge . coefficients CDF 
and CDX *) . , 


Fuel Injection AP 
LOX Injection AP 
C*-Eff iciency ' 

^sp 

Combustion Chamber Pressure 
(injector face, static) 


The listing of the deck setup for input data is shown in Table 4-1. 


*) For an engine for which analytical or test results are available, 
APFUEL and APOX can be used to calculate the discharge coefficient 
for the orifices. Using Equations ( 7-6? ) ( 7 - 70 )ft*om Section 7.3, 

Part 2: 


WINJF 

CDF = 

ttDF^ I^P 

4 

and 


WINJX 
CDX 

NOX 

4 


/2gRH0X APOX 


2gRH0F APFUEL 


For definition of symbols see Sections 2.3 and 2,4 
(g ;= gravitational acceleration) . 


• i 
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AFl = 150.* 

AFII = ?00.» .^.. . 

AFIIl = .00069V 
triv = .03^St ' 

AFV = 40. 5» 

AFVI « 3S.* 

AFVn = 66,* 

BFI = ).0* 

9FII = O.Ot 
BFIir- = l.A* 

BFIV ^•0? . - 

PFV * 1 .B6t 
PFVI = ?*0» 
prvn = i.i6f 
crn s .It 
Am = o.t • - 
AF j?=o.» 

AF12 = l.t 
Ari3=0.# 

AF13 = l.t 
AFlArO.t 
AF14 = .Olt 
AFlS=0.t 
AF15 =. .Olt 

AF20=0.t 

AF?0 = l.t 

AF21r0.t 

AF21 = l.t 

AF2?=0. t 

AF22 = l.t 

AF23=0.t 

AF24 = O.t 

AF?5=0.t 

AF^^F- = l.t 

AF>6=0.t ‘ 

: AF2A = l.t . 

' AF27 r O.t 
AF27 = l.t; 
AF2P=0.t 
AF?H = l.t 
AF?Q=0.t 
AF29 “ l.t 
AF31=0. t 
AF31 = 2.t 
AF3?=0.t 


B 


TABLE 4- 



Ctiiiiiiq 

1 i 

f '■ 

1 ENGINE Tt'PE NO.l 1 

m 



ITEST CASE 

^ TH.^UST = 15000 POUND/ 
* 02 / H2 propellant ‘ 

m 

AF4J = 

.009. 

m 

AF42=0 


m 

m 

AF42 = 

• 12S. 

1 optimization 9UN 1 

m 

m 

AF43=0 

• . 

J — 1- 

m 

AF43 = 

.125. 


AF44=0.* 

AF44 o S.6?5t 
AF4f= 0. t 

AF45 = 4.16t 

AF46=0. t . f 

AFSl =0. * 

AFSl = l.t ■ 
AF5?=0.t 
AFF2 = l.t 
AF4l=0.. 

AF61 = l.t.- 
AFh?=0.t 

AF^.? =1.. 

AF7l=0. 

AF71 = 4.0t 
AF72=0.t ' - 

AF7? e 2,0t 
BF7l=0.t 
BF71 = 5.06t 

*cQ»^npv 

OF * .l?92. 

DX s »0Fi4, 

NF = ? 1 6 • t 
K'OX = 216. t 
VOLF = lO.t 
VOLX = lO.t 
LF s .06* 

LOX = .40t 
LC = ll.??t 
OC = 10.?bt 
CF ss 5.06t 

s - 

SiCOVMlS 

th;<ust .= isooo.t..; 
WT = 33.7?t 
i/FCOOL =' O.Ot 
WFFPF.B .= O.Ot :: 
WXPPEB = O.Ct 
LI = A.Ot 
OTHPTI e 5.14. 
DFVFL = .1292. 
oxvEu * .not 

ANGlMP=0.n. 
CnFs.6'3?t 
CnXr.6C5t 
Lfi?0Xt:7.7S. 
LF.50Fs2.75. 
RF1=0.0. 

«?F?=0.0t 

LP s 0 . . 

CPFcl.OTt 
CPXs.410. 

RHOF « .OOOlSlt 
TLOF s 349.. 

PMOX s .0404. 
STGLXs. 0000453. 
MUX K .00001 t 
TLOX^a 252. t 

m 
m 
m 

■ 

t 

. Complete Listing of the Deck Set Up 

for Input Data of Sample Engine (continued) 
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Atiiiiiiiia 


t 

KOLFX=3^.0* 

LA'VOK = 3‘^.0» 
TC^x-?7fl.6. 

CHIF = 

fLiNH = .Oftf 
ruA»EA = 2.7A, 

CHlX = .OA, 

XtWZ = .AO* 

XLAt?£A = 1,19, 

ETYP =1., 

PTYP =!.♦ 

ICOax = 1* 

INJC^^P = ?f 
lEOPrl. 

K0x = u 

istatf=i, 

istatx^?* 

ISTOPFrO* 

1ST0PX=^0* 

PCPC(l)sA*C-.Of 

PCflF(l)=A<^L‘.0> 

PCRX(l)=4<»O.0t 


KUF=0.1S9£-6, 

LAt*P.PF = 19?.0, 

SK>LF = O.Ol , 

TC?F=SO,5* 

LA^1F12=192. 

F23TBS=1.9, 

NTAi;ni = .a 5 * 

NTAUi'<3=.7f 

$ 

V^LPODT 

I=>f<00T = .TR'JF. * 

NJ?F20 = l*NTF?C;s3*NRF2l=3»NTF21«3f 
NPF?6=S*NTF26=5, 

N'>F?6ri,NTF?h = 3» 

N9F?#>=2,NTF2rf = 2» 

NF?7=7. KB?8=S, 

N^IF ?7 = 3 *NTF 27 = 3 i 
ROOTS=.TkUE . , 

ROOTS=. false., 

IG^FPO-. false., 

1GPF21=. false., 

jnr/F?A=. FALSE. , 
lr,PF?7=. FALSE. , 

IG»F2A=. false., 

H'pr?n=.T?nE.. 

2l = .TPl»F., 
lG?.7PS=.TRUr.* 

IGPF?7=.TPUF., 

ir>-^F?>i=.TPUE., 

PCU7PS-I0.C, 

RCUT27=1C.C. 

AI.PKi?q=10.0. ALPN27=lO0.O, 

ALRI>?*» = 10.0, ."•LRN?7 = 1S.0, 

&oF20Pn .37S?,-n .3752* 

GRF20I ( 1 ) =0. *“n.3;.A2,13.:5>-A2* 

GPF21PM -.^OOOSS. 

GRF21 I (1 )=1 .59S77,«1 .5vS77 , -3. 9RCA 1 2E-09, 

GRF?f.Cf ( 1 ) S-. 3 AS 095 . - . 21 F l.‘jA ,.".21 ol ai H* j 5227 . 

6RF2M ( 1 ) =9,h99e75E-0B,P.si9a31 , -2.99^.32 , -3. tt3S8 1 2E- 07^ 


i T 

- GRF27I(l,l)=l.S<‘129, 

S riRF?7R{l*2)=0., 

- GWF?7I(l*2)=-l.fAi29. 

~ G9F27R(2,1)=.03A99, 

Z G9F?71<2*1)=3.0535, 

Z GRF27R (?,?1 =.03A977, 

Z GRF271 (2.2)=-“3.C535, 

; GPr?7R(3*l >=0.* 

Z GRF271 (3» 1 ) =3.C32, ' 

- GMF27R(3*2)=0.“* 

• GRF27I (3,2) =-3.‘h32, 

- GPF27R(4,1)=.0)A,' 

1 GPF27I(A,1)=7.016, 

3 GRF?7R(A,2)=.0l*', 

■ GPF?71 (4,2) =-7.C16, ' J 
Z G«F?7R(?,1)=-. 013549, 

Z GOF271 (S,l)=5.2934, 

Z GRF27R(S*?)=-.Cfl3&^6♦ 

z GPF27I (5*2) =-S. 2934, 

3 GMF?7fl(ft,n=0-, 

3 GRF?7I (G, 1 ) =><.526, 

3 GRF27R(6.2)=0., 

- GRr ?7I (6,2) =“8.S25, 

Z GRF27R(7.1 ) =-.?C95, 

: GRF27I (7, 1 ) =6.RL184, 

3 GRF?7R (7,2)=-.?C'9S, 

3 GRF27I (7,2) =*6.F01fi4, 

3 GRF?P,I (1,1 )=.l3rSJ4, 

- GWF28R(1 ,1) =1 .76662.,' . 

• GRFPPI (1 ,2) =. 1 36B24, ’ 
3 GRF?HR(1 ,2) =-l .7B661V 

2 GRF?Pi{2.n=.27c302, ■ 

3 GRF2FP(?,1)=3-?p=^39, • 

3 GRF28n2,2) =.27426'4,- 
3 GRF2»^P(?.2) =-3.29«.36r 
3 GRF2HI(3,l) = .273SC6i'i 

• GRF?PR(3,D=4.3F735,;' 

Z 6RF2«J(3,2)=-273462, 

2 6RF28R (3,2) =“4.3B732:, . 
2 GWF2P I {<*,1 ) =.275406,/ 

2 GRF23R (4, 1) =4.03222,- - 
2 GRF281 (4,2) =.275363,, 

2 GRF?ftR(4, 2) =-4.03220, 

3 GRF2/3I (5,1 ),= .C?6 449.9, 

• GRF?6R (5, 1) =.027628, 

- GRF2BI (<,2)=.02H4498, 

Z GRF2PR(S, 2)=-. 027623, 

Z S 

3 SCO'^TPw 
3 IRwU = l. 

3 IRWU=0. 

- !Rwi2=i, 

- IRv/12=0, 

Z IRW13=1, 

2 TRwn=o, 

3 iRwi 4 =i, 

3 IRv;14=0. 

3 TPW15=1, 

• 1PS*'1S=0, 

2 IR*^?C = l , 

2 iPv*?o=o, 

3 IRt.*21 = l, 

j TRw?l=0, 

• 1RW22=1, 

2 1RW??=0, 


▼ 


TABLE 4-1 . (continued) 
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IPW?^30. 

IRWPS=1 
19v 2S=0. 
lPiC36 = ?. 

1RW?A=0 
IR'^PTs?* 

IP--^7rl 

IPw?7=0 

lP/:>rt = l 
IPVP^* = 0 
IR-4?<;=1 

IPV?«?=0 
IPwil=l 
lPw3l=n, 

IPW V = 1 
IPw3?=n, 

IQr.-41rl 

IPw<»i=o 

IPw/i3=i 
IPW4?=0 
IPWA3=1 
JKV4.3=0 

IPV‘46 = 0 
IPW<.Srl 
IPV45=0 
IPW4^=1 

lpV:.6=0 

1RWS1=1 
JRwSl=0 
IR>’M = J 
1P>.M=0 
1R'.%'71 = 1, 

I«W71=0» 

IPv;viPsl* 

lPwVftP=C» 

TRh’Pal = U 
IPfr'PALrOt 
JPwt50*lt 
IP>*L51J = 0. 

IPPOOT^.TPuE.t 
IPMISC*.TRUE.»' 
$ ‘ - ; 
scpvni 


ipiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii| 


ALPHAn>=.l?ni*.071?l'i4,??>9.?57,?AV.:>S7*U.5i<ri3a*3 0.'^S73 
ALPHA(7)=.0(S653b^».37S51f l.C0»9.rt707S,7,o<»£y71t«t.957c,«»* 
MAXJJJ =0* 

S 

OOCOOOOOOOOCOOOOOOOOOO 

OOOOOOOOOCOOOOOODCOOGO 


V 


\ 


\ 


\ 


\ 


\ 


IPMAf/L = ,F4LSE.» 

S 

y Tftc‘<;op 

IWFPf.T=:20« 
tUPI^=?• ITFIN=C» 

MUmopt^ 2 , KETH 0 P(l)= 9 * 2 f 
►•AXPPSsSf • 

IPr,AlVaU 

I«<iPCH3U IPACYCe;, 1PFCYC»1* 

MAXJJJ = 2f 

«AXJJ.I=:30% 

MaxjJJaSS. 

PAXJJ.J=l00t 

ALPLO<l) = .»'57%.D^3tl62.,lf>?.%7*5»7.Si.0AS*.3f4.*ft.42»7.«il.3.79* 
ALPHin>s.l6J5*.l05.^70.f27C.*J?.S*'12.5f075t.S») **14.04, 12.86* 5.0599^ 


biiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii 


, (End) 
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TABLE 4-1 



4.2 Data Output for Test Case 


Computed running conditions for this engine from the AUTOCOM code were 


Fuel Injection AP 

82.9 psi 


LOX Injection AP 

48.5 psi 


C* Efficiency 

see Table 

4-3 

Isp , 

444.8 Ibf 

sec/ Ibm 

Chamber pressure (injector 
face, static) 

372 . 8 psi 



It is assumed that the low chamber pressure computed results from the 
ratio of specific heats employed for the propellant combination (Y = 1.2505) 
and the combustion ten^erature = 5722“R). A more complete set of 
computed engine running conditions is presented in Table 4-2. 

4.2.1 Specific Combustor Characteristics ; Specific combustor charac- 
teristics for the nominal engine are presented in Table 4-3. The average 
engine characteristics resulting from the selected, specific combustor charac- 
teristics are shown in Table 4-4. The rating value resulting from the 
selected average engine characteristic weighting fa:Ctors is also presented 

in Table 4.4.' 

It should be noted that the specific combustor stability characteristic 
consumes almost all the computational time required for the. evaluation of 
a combustor. This is' due to the time consuming complex characteristic 
equation solutions required for chugging (F20 and F21) , longitudinal time 
lag stability analysis (F26) , transverse time lag stability analysis (F27) 
and the Lewis response function stability analysis (F28) . Table 4-5 presents 
a summary of the characteristic stability equation roots for the nominal 
engine. The least stable root is obtained from the transverse time lag 
analysis (F27) using a value of = 3.0543. 

4.2.2 Stability Roots . Some difficulty was initially experienced in 
computing the nominal engine stability characteristics for the time lag 
analyses. The AUTOCOM program assumes a value of the Reardon interaction 
index, n, of 0.5 for the longitudinal time lag analyses and 1.0 for the 
transverse time lag analysis. With these interaction index values, the 
nominal engine was found to be slightly unstable in two of the transverse 
modes. Table 4-6. A sensitivity study on the effect of interaction index 
value was undertaken. As a result, an interaction index value of 0.90 was 
subsequently utilized in all transverse time lag analyses and an inter- 
action index value of 0.45 was used for the. Ion gitudinaL; analysis . These 
values were used to obtain the time lag analysis roots shown in Table 4-5. 
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ND.l 

Tr.bT CfcSf 

Tri^UST = l*iUOO P0U>10 
0<i / M? P><OPtLLAUf 
cpriMi^Anow ^UN 

OENi^^AL F.NOINc P#»3 amF.1^^S 


COM^USTOP ThPUST FOf^CE. 



THRUST 

= 

15000.0 

PPnPFLLflNT SPtCIFIC IMPULSE 



ISP 


AAA, 6^0 

CHt'«REP PRESSURE AT INJECTOR HEAD 



PCI 


37?. 765 

total pl'opellant flow rate 



WT 


. .33.7200 

FUEL WEIGHT FLOW t?ATE 



WF 

= 

5.S64A5 

OXIOI7.ER weight FLOW RATE 



WOX 

= 

28.1556 

COMBUSTION TEMPERATURE IN CHAMBER 



TCOMB 

= 

576A.91 

RATIO OF SPECIFIC HFAT OF COMBUSTION GAS 


SPHEAT 

= 

1.2A9AA 

•..GAS CONSTANT OF COMBUSTION GAS (ft, lb, 

./Ib.OR) 

RGAS 

= 

127.127 . 

ioeal thrust coefficient 



CFIDEL 


1 .93928 ‘ 

ACOUSTICAL LENGTH OF CHAMRFR 



LOUM 

= 

8.81 57A * 

KOLECULtp wEIGhT'OF COMBUSTION GAS 



MOLWT 


12.1A53 

MEAN RESIDENCE TIME OF GAS IN ChahrER 



THETAG 

= 

8.567267E-I 

SPEED OF SOUND IN CHAMBERS 



CS 


65133.9 

... COMBUSTION chamber »^aCh NU^'BER 



MC 

= 

;20t068 

INJECTOR PRESSURE DROP FOR FUEL 



OELPF 

= 

82.90A9 ' 

’ INJECTOO PWESSUPE DROP FOR OXIDIZER 



DELPX 

= 

A8.A525 

COMBUSTION chamber VOLUME 



VC 


681.565 

AVERAGE VELOCITY OF GASES IN CHAMBER 



VELC 

= 

13096.A 

FUEL INJECTION VELOCITY 


- 

vruEL 

= 

1301 3, or 

OXIDIZER INJECTION VELOCITY 



VOXIO 


339.511 

TOTAl A«>EA OF FUEL INJECTOR ORIFICES 



TAF 

= 

2. 831 BA 

TOTAL AREA OF OXIDIZER INJECTOR ORIFICES 


TAX 


1.19702 

CHAVREM length TO VAPORIZE 50 PER-CENT 

OF 

fuel 

L5FUEL*= 

0. 

.. chamber length to VAPORIZE 50 PER-CENT 

OF 

OXIOIZER 

L50XIO 


.636932 


* Fuel in a gaseous state 


TABLE 4-2. Nominal Engine Running Conditions 

(Dimensions are in in., lb., sec, °R) 
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NIIMPFW 3)?> 


Ttpe »jo*i 

itvr 

T*^yi‘ST = lbo30 PO.Pvj) 
0? / M? P^^O^^ELLANT 

Optimization »^un 
sprcinc COf^BUSTOP Cda-^ ACT£;»/IST ICS 
PE*? CENT MASS FUEL VAPOfJIvtr) 

PF^ CENT mass of OXIOIZEk VAPO^TZEO 
C» FFFTCIENry MIXING MOPEL 
C* FFFICIEnCY PULSFO COMhOSTOPS 
C* FFFICIENCY MOM-PULSF.D COMmUSTOPS 
FUEL system CHUGGING DECAY RATE 
07 !IDI 7 ER SYSTEM CHUGGING DECAY RATE 
PULSED INSTABILITY CHAPACTFP I ST IC 
NON-PULSED INSTABILITY CHAPACTEP I ST I C 
OY?<EHA FUEL STABILITY DECAY KATE 
DYKEHA OXIDIZER STABILITY DECAY RATE 
stability longitudinal TIME LAG 

stability transverse tipe lag 

S I AfiiL 1 I I L*^c KC SkunSl f' ONC 1 IOn 

stability priem linear analysis 
FUEL PPESSUOE OPOP ChAPACTFR 1 ST I C 
OXIDIZER PRESSURE DROP CHARACTERISTIC 
FUEL PLUS OXIDIZER HOLES CHAKATER I ST IC 
OXIDIZER DOME VOLUME CHARACTERISTIC 
FUEL DOME VOLUME CHARACTERISTIC 
OxlOIZER HOLE LENGTH CHARACTERISTIC 
FUEL HOLE length CHARACTERISTIC 
INJECTOR TYPE COMPLEXITY CHARACTERISTIC 
INJECTOR LFNGTH CHARACTERISTIC 
CHAPHER DIAMETER CHARACTERISTIC 
MIXTURE RATIO CHARACTERISTIC 


Fll 

100. (Gaseous ) 

Flp 

= 99.4122 

F 13 

= 100.000 

F 14 

= 91.3821 

FIS 

= 73.7665 

F 20 

= -1036. 

F 21 

» - 261,5 

F ?2 

e - 37.4819 

F 23 

s Not Computed 

F 2 A 

a Not Computed 

F 25 

» - 2964.82 

F 26 

« - 1681.31 

F 27 

« - 412,843 

^28 

a - 360.322 

F 29 

c - 46427.6 

F 31 

= 82.9049 

F 32 

a 48.8525 

FAl 

= 431.000 

FA 2 

s 10.0000 

F 43 

a 10.0000 

F 44 

a .400000 

F 45 

= 6 . 000000 E -02 

F 46 

a Not Computed 

F 51 

a 2.18268 

F 61 

a 2.00000 

F 71 

a 5.05990 


TABLE 4-3. Nominal Engine Specific Combustor 
Characteristics 
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316 


ENGINE TYPE NO* I 
TEST CASE 

thrust = 15000 POUND 
02 / H2 propellant 
OPTIMIZATION HUN 


AVERAGE ENGINE CHARACTERISTICS • 

performance characteristic 

FI 

s 

.536361 

RATING component 
85.4541 

STABILITY CHARACTEHISTIC 

FII 

s 

-37.4519 

4.71208 

PRESSURE OHOP characteristic 

Fill 

s 

71.4208 

1.45876 

INJECTOR COMPLEXITY CHARACTERISTIC 

Fiv 

= 

8.37860 

3.11377 

LENGTH CHARACTERISTIC 

FV 


2.182M8 

173.002 

Chamber diameter characteristic 

FVl 

s 

1. 00000 

35.0000 

chamber mixture patio characteristic 

FVII 

s 

4.000000E-08 

1.729552E-07 


COMBUSTOR RATING s 2'i7.780 


*O00tf0»«00««« 00000 »« 


TABLE 4-4. Nominal Engine Average Characteristics 
and Rating 
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NUMBER 316 ENGINE TYPE NO.l 

TEST CASE 

THRUST s 15000 POUND 
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TABLE 4-5. Stability Roots in Time Lag Analysis 
C F26: n = 0.45 ) 

C F27: n = 0.90 ) 



A second point should be noted regarding the stability roots. No root 
is found corresponding to the third value of Svh = 3.8317 in the transverse 
time lag analysis. The missing root can be found by varying the initial 
guess value in the complex plane for this particular root. Following this 
procedure the missing root was found to be at the point (-.59424 ± j3.5144) 
where the Reardon interaction index was 0.9. The root is thus highly 
damped. 

When commencing an optimization study, side analyses of the above type 
may be required to locate particularly difficult roots. This procedure 
should also be followed whenever the root imaginary part is not approximately 
equal to the corresponding value of Svh transverse time lag analysis 

and the Bessel argument, m, in the Lewis response function analysis. This 
point is discussed further in Section 4.3.1 s 


4.3 Optimization Computations for Test Case 

Optimization computations on the sample engine were initially undertaken 
using all specific combustor characteristics and all stability roots. 

However, it was noted that the design variable perturbations introduced little 
change in the computer time consuming stability equation roots. Accordingly 
the combustor analysis was divided into two classes of computation. These 
were an approximate analysis which considered fewer (possibly none) of the 
stability roots and a oorrplete analysis in which all stability roots 
were computed. It is emphasized that the approximate analysis is only 
approximate in that the calculation of the less significant stability roots 
is omitted. Clearly, by a judicious mix of complete and approximate analyses 
the total elapsed computer time required for the definition of an optimum 
engine design can be drastically reduced. 

4.3.1 The First Twenty Iterations . Following initial experimentation 
using all stability roots, the engine was subjected to twenty design iter- 
ations using all specific combustor characteristics. An approximate analysis 
mode was employed which considered only the relatively rapid calculation for 
the longitudinal time lag analysis (F26) and the transverse time lag 
analysis (F27) for the single Sv^ value of 3.0543- (the least stable transverse 
time lag root). This approximate analysis permits both longitudinal and 
transverse stability characteristics to be monitored. 

Initial and final stability roots from this optimization calculation 
are presented in Table 4-7. It can be seen that little change has occurred 
in the stability roots. The trend is to increased stability in the less 
stable transverse mode and to less stability in the more stable longitudinal 
mode. It may be noted that the transverse time lag analysis of Table 4-7 
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' A. SELECTED STABILITY ROOTS, NOMINAL 

•>uHa£0*3lfi TY^t NO.l 

TtST case 

THajiT s ISOOO POUND 
02 / PPOkEllANT 
0?TI:4I2AT10N PUN 


STaRILITY poot su^kasy 

ROOTS FRoV SECr?-i (STAnlLlTY LOnGITUOInAl’T IM£ LAG) 
COAXIAL injection 


ROOTS SECF27 iSTArlLlTY TRANSVERSE TIhE LAC) 

COAXIAL injection IFRcO ■ 1 SNUH »• 3,05A30 


ENGINE 


IntERhEO'IATE combustor output jjj 


DECAY RATE 
-.325761 
-.227565 
-.22756A 

DECAY Rate 
-3.257S76E-02 
-1,22075 


FREOoEk’CY 

S.5a0n34£-03 

-3.01A«:2 

FREQUENCY 

3.AS2U 

-1.52265 


I 


i , • ' ' 

B. SELECTED STABILITY ROOTS AFTER TWENTY DESIGN PERTURBATIONS 

31<* F.NUTNE TY^E nO.I 

TK 5 T CASE 

T"iyi;^T = 15000 POUND 
02 / »^‘* 0 --tLLANT 

OPT I M /A r ION PU.N FINAL COMBUSTOR OUTPUT JJJ=21 


STAt'lLTTY -001 SUAMAKY 

ROOTS SrC^pA <ST'*vli I TV LONG I Tui) I N AL TIm£ LAG) 

COAXIAL injection 


HOJTS F-Yim SrCF 27 «ST-*rIl ITY TRANSVERSE TIkE LAG) 

CoA<I^l. INjtCTjt.'* I^RtJ ■ I 



- 

DECAY RfiTE 
-.315557 
-.223750 
-.2237A8 

FREOflF'-tCY 
4 .2 ^Rs>A 1?-07 
3,r,iM<) 
-3.01520 

SNUM ■ 

3.05430 

OtCAY RATE 
-A.67317AE-02 
-1 .2A530 

FREOilcMCY 
3.f»6709 
-1 .57179 


TABLE 4-7. Selected StabiHty Roots for Nominal Engine,. 

and Engine after Twenty Design Perturbations 
( F26: n = 0.45 ) 

( F27: n = 0.90 ) 


31?^ ENG INK type; NO.1 

TF.ST CASE 

THinisT = isono pound 
u? / h2 propellant 

, 0Pri«I2ATT0rj RUN 
average engine CHAPaCTEwISTICS 
PF»'FOP*-tANCc Chap act F» I FT I C 
STA'^IUTY Ch as^.CTF.t^ ISTIC 
PJESSIJ^F PPOp ChaPACTEPISTIC 
JN.'KCTOt^ C04PL«:XITY CHAPACTKPISTIC 
LF'^M h CHA'-'ACTr i'lST IC 

f-IAf^ETFk CMAKACTEPISriC 
ChAkPcP v1y7UP£ ratio CKA-aCTEPISTIC 


Rating component 


FI 

= 

.344513 

51.6770 

FIl 

s 

-S6&.349 

4.1491905' 

Fill 

= 

75.3314 

1.64971 

FIV 

= 

«.902b4 

3.13079 

FV 


2.20692 

176.563 


3 

1.0H42<J 

41 .1487 • 

FVII 


7.C9733BE-05 

1.01S815E 


« O »««> 4^ » l» « O' « 4» «-«»«■««« « 


C0MPUS70R RATIfxiG = 274.170 


TABLE 4-8. Rating after Twenty Design Perturbations 
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considers two solutions to the stability equation with Svh 3.0543. These 

are the true solution with the frequency approximating Svh> ^ spurious 
solution with the frequency approximating tt/ 2. These spurious solutions 
with a frequency approximating tt/2 are often encountered in the time lag 
analysis. If the true solution is not obtained on the nominal engine evaluation 
and the spurious solution is obtained, the AUTOCOM program will "track" the 
spurious root. Hence, the analyst must take care to insure that the correct 
roots are found on the nominal design before embarking on an optimization 
run. This point is also discussed in Section 4.2.2. 

The engine rating after twenty design perturbations and the corres- 
ponding average engine characteristics are presented in Table 4-8. It 
can be seen that based on the selected average engine characteristic weights 
which provide the rating in the form of payload lost, a gain of 23-1/2 
pounds of payload has resulted when compared to the nominal design of 
Table 4-4. It can also be seen that the average stability characteristic 
contribution to the rating is now negligible and that payload is being 
gained primarily by reduction of the performance characteristic penalty. 

Pursuing this payload improvement, it can be seen. Table 4-9, that the perfor- 
mance improvement stems from F12, per cent mass of fuel vaporized, and from 
slight improvement in C* efficiencies for both pulsed and non-pulsed combustors. 

4.3.2 The First Hundred Iterations . Following the first twenty design 
iterations discussed in Section 4.3.1, the optimization problem was restarted 
without any stability analysis; and 100 successive design perturbations were 
introduced. A combination of the uniform random ray and pattern searches 
were employed. Reference 4-1. The approximate analysis employed completely 
neglects the stability characteristic. The rationale for this approach was 
the negligible stability characteristic contribution to the engine rating. 

Table 4-8. This table indicates that the stability characteristic affects 
the rating in the twenty-fourth significant figure. This is well below the 
accuracy of the CDC 6600 computer which, with sixty bits, is able to provide 
approximately ten significant decimal figures. 

The nominal engine rating without the penalty of all stability charac- 
teristics (4.7 pounds. Table 4-4) is 293.1 pounds. After 100 successive 
design perturbations introduced through the References 4-1 and 4-2 multi- 
variable search program, AESOP, the rating is reduced to 210.7 pounds. 

Rating convergence is illustrated in Figure 4-1. Convergence behavior of 
the combustor design variables is illustrated in Figures 4-2 through 4-4. 

The combustor design variables were allowed to fluctuate by plus or minus 
twenty-five per cent of the nominal values in this study. Two of the design 
variables, the chamber diameter and the niomber of fuel orifices (which 
equals the number of oxidizer orifices) are practically on the lower and 
upper bounds permitted. 
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:{jV«rP. *115 \ TY?c wo.l 

IRSl CASK • 

' TH.^UST = li>O00 POUND 

up / h? KPvU^KLLANT 

. • . ‘ ^ OPT IMI?**! ION kUN 

S=*'CI*^1C CO- huSTOP CMAyACrePIST ICS 
P F P Cr Vi 7 ■ i 5 s r i/r.L V A P 0 i 7 F! D 
PFP Cr'Jl MASS OF OXIOI’Fr/ VoPOPIZF.O 
C.*-rFFlClENCY P.OLEl ; 

C* EFFICIENCY PV?LSF.O CONMwSTOHS ^ 

C« iTFFICIENCY NON-'puLSED C(»f-*MUSTUPS 
FUEl. SYSTEM CrtUOOING DECAY PATE 
OZlr>I7EP SYSTEM CntiGGlNG r*EC6Y kaTE 
PULSF.n I^:STA9IL1TY CHAPaCTEkISTIC 
NON’-, PUL Sc 0 INSTABILITY CmapaCTlW IS T 1 C 
DYKEMA FUEL STABILITY DECAY PATE 
•0YKE''*A OXirUZtP. STaJ^ILITY" decay KATE 
"STAHILITY iONOiniDINAL TI«E LAG 
STABILITY TBAUSVEkSE TIME LAO 
cT^oTt jTv r i>r b:,c j/^.vc c 

STARILI TY ,PBjt -4. t iNEAft ANALYSIS ; 

FUEL PPESS'Ji?E CROP ChakaC f r K 1 ST I C 
OXini7EP PRESSURE DROP C*^AP aC TFP i ST I C 
FUEL PLUS OXIUI7EP HOLES CH APACTE» I ST IC 
OXiniZEP po’E VOLUME CHartACTc.PlSTlC 
.FUEL -DOMf VOLUME ChaKaCTEk I ST I C . 

OxirU7FP HOLS. Lb'ir.TH CHARACTERISTIC 
FUEL HOLE Length ChaK ACT E^ I ST I C " 

INJECTOR Type COMPLEXITY CH akACTER 1ST I C 
INJECTOR LENGTH Cn AP«C T F R 1 ST I C 
charrep ima'^'Eter Character 1 STIC 

MIXTURE PAT. TO CHAMACTEKISTiC . . j 


r 



Fll 

S 

100. 

F12 


1 00.000 

F13 


100.000 

F14 


91.7069 

F15 

S 

73.8398 

F20 

= 

0.' * 

F21 


0. 

F22 

*; 

0. 

F23 

■ = 

0. 

F24 

= 

0. 

F2S 


0. 

F26 

= 

-1636.16 

F27 

= 

-566.349 

F2S 

r 

A ^ 

F29 

s 

0- 


= • 87.9A79 
= 50-0986 - 

= 433-943 

* 10.02i?3 

s . 9.97392 
= .399630 . 

= 6.005603E-02 

s Not Completed 

= 2.20692 

« 2.08429 

= 5.05579 


F31 

F32 

F41 

F42 

F43 

F44 

Fas 

F46 

FSl 

F.61 

F71 





TABLE 4-9. 


Specific Combustor Characteristics after 
Twenty Design Perturbations 
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Ojily fgg and f2'j stability 
roots were computed 
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Figure 4-1, Combustor Rating Convergence 
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Figure 4-3. Figure 4-4. 

Combustor Design Variables Convergence Combustor Design Variables Convergence 


The final rating and the characteristic components to the rating are 
presented in Table 4-10. Final design variable values together with the 
search limits employed are tabulated in Table 4-11, From Tables 4-4 and 
4-10 the rating changes associated with each characteristic are seen to be 


Performance Characteristic 
Stability Characteristic 
Pressure Drop Characteristic 
Injector Complexity Characteristic 
Length Characteristic 
Chamber Diameter Characteristic 
Chamber Mixture Ratio Characteristic 

Total gain; 


20.9 lbs., gain 
Not considered 
0.14 lbs. , gain 
0.72 lbs., loss 

36.7 lbs., gain 

26.7 lbs., gain 
1.32 lbs . , loss 
82.4 lbs. , ... 


The total rating gain of 82.4 lbs. produced by the optimization process 
of the AUTOCOM code ignores any stability characteristic effect. To assess 
this effect, a con 5 )lete analysis was performed using the Table 4-11 vector 
of combustor design variables. The rating resulting from this complete 
analysis is presented in Table 4-12. The associated specific combustor 
characteristics are presented in Table 4-13. The stability characteristic 
produces a rating component of .16 pounds, a 4.55 pound in 5 )rovement over 
the nominal engine stability characteristic*. Comparing the final rating 
of 210.84 pounds. Table 4-12, with the complete nominal engine rating of 
297.78 pounds. Table 4-4, the total rating gain obtained in 100 design pertur- 
bations is 86.94 pounds. It is interesting to note that despite the use of 
an approximate analysis which resulted in the stability characteristic being 
ignored, this characteristic nonetheless improved during the 100 design 
iterations. Elapsed computer time for the 100 iterations, the final complete 
analysis, and the initial complete analysis was 250 seconds on the CDC 6600 
computer. 

4.3.3 A Note on Stability Roots After 100 Iterations . The complete 
stability root set obtained after 100 iterations is presented in Table 
4-14. It can be seen that the second frequency corresponding to S^h “ 3.0543 
is missing. This root was the least stable on the nominal engine, -Table 4-5, 
but became more stable in the first'' 20 iterations of Section 4.3.1, Table 
4-7. Accordingly, a search for this root was initiated to confirm the 
stability improvement over 100 iterations. The root was located as a non- 
conjugate pair at the points 

= (-.443094 + j2. 73775) 

and 

Z2 = (-.552853 + j3. 61289) 
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31^ £Nr,I\£ NO. I 

f I T£s: CAS£ 

Th.»US7 ■ 1^000 POUND 
0 ^ / P^^CPtLLANT ’ 

- optimization Sun 
AVEPiGE Engine characteristics 


prRropyANCE CHAPAcre»isTrc 

FI 

m 

stability characteristic 

FlI 

- 

PPESSUPE Dm'OP CHAKACTckISTIC 

Fin 

o 

INJECTOR complexity ChAPACTER ISTIC 

FIV 

o 

LENGTH CHARACTERISTIC 

FV 

- 

CHAMBER DIAMETER CHARACTERISTIC 

FVI 

tt 

CHAMBER MIXTURE RATIO CHARACTERISTIC 

FVH 

m 


C 

NOMINAL RATING 



rating co».ponent 

COJJlPONSNT, TABLE E 

.397158 

59.5737 

80.U5l*l 

0. 

OO.QOO 

u . 71208 

67.5013 

i 

1 .3S3Q8 

I.U 9876 

9.856fl0 

3.63768 

3.11377 

1.9203B 

136,319 

173.002 

•4B6520 

6.2S4S7 

35.0000 

3.426861E-02 

1.31833 

1.729552 X 10 ' 

COMHUSTOR rating 

* 210.6R7 



NOMINAL COMBUSTOR RATING » 297.760 


TABLE 4-10. Rating after 100 Perturbations, 
Approximate Analysis 


• t 


M»HRER 219 £Nf‘*TNE TYPE NO.l 

■TEST. CASE 

. THRUST *> ISnOO- POUND 
02 / M2 PROPELLaJ^T 


optimization run 

final values of tmk design variables 



{ O^^High 

oC Low 

Noninal 

diameter of fuel orifices 

OF 


.IPlllO 

j .1615 ■ 

.097 

.129 

diameter or OXIDIZER ORFICES 

-OX 


7.121544E-02, 

.105 

.063 

,.08U 

NIJHMEP of fuel orifices 

NF 


260.257 

|, 270. 

162. 

216. 

number of oxidizer ORFICES 

- NOX 


260.257 

270 .. 

1 

162. , 

216. 

VOLUME or fuel manifold 

VOLF 


10.5458 

1 12.5 

7.5 

10. 1 

volume of oxidizer manifold 

VOLX 


10.3573 

1 12.5 

7.5 

10. 1 

LFNGTH of fuel ORIFICES 

LF 


6.653641E-02 

1 .075 

. 01*5 

.06 1 

length of oxidizer ORFICES 

LOX 


.37MS10 

1 - 5 . • 

.3 

.u . 1 

element Type 

etyp 


1 .00000 

1 

1 . 

1 

length of CmAmPER 

LC 


9.A7075 

lU.ou 

8 . 1*2 

11.22 1 

Cha.mber*OU.mETEH 

DC 

3 

7.64071 

12.86 

7.51 

10.28 j 

MIXTURE RATIO tOXIOIZEP/FUtL) 

OF 

3 

4.06744 

1 5.0599 

3.79 

5.06 1 

FINAL COMbUSTOR RATING 

8 

210.607 





KOMIHAL COraUSTOH lUCEG >= 297-780 


TABLE 4-11. Design Variable Values after 
100 Perturbations 
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3? a 


TY'»t NO. I 
TcST C*SP 

ThHJST ■ ISOOO POUND 
02* / P*<OkKLLANT 
0»T1HI2ATI0N PUN 

AV£C.AOt ENGiKr. CHABaCTEP I ST ICS 

PE»FO»‘**nCE ChapaCTcPISTIC 

STAqii.ITY ChAhaCTcPISTIC 

ppr<;<;jor OPOP C^A’^'-CTEBISTIC 

INJFCTOO C0«PL^J«ITY CriaPACTEPISTIC 

length CmapaCTFPISTIC 

CHAPPEB DIAvtTfP Cha-«aCTF.PIST1C 

CHAMREO MIXTURE PATIO CMAiiACTEPISTlC 


KOVH?AL RATING 

rating cosponent coMPosisrr, table hi 


FI 


.397lSft 

SP.S737 

80.1*51*1 

FII 


-71.AIH3 

.ISnPAU 

1*. 71208 

Fill 


A7.S01 1 

1.3S397 

I.U 9876 

FIV 


9.e56H0 

3.e.37AH 

3.11377 

FV 


\ .V?03^ 

136.319 

173.002 

FVI 


.AbSSIQ 

A.?R454 

35.0000 

FVII 


3.a2'69a1E>0? 

l.:slH37 

1.729552 



COMhUSTOR mating • 

210. R45 



NOMINAL COMBUSTOR RATING « 297.780 


TABLE 4-12. Rating after 100 Perturbations, 
Complete Analysis 


MJMdER 3?0 tSftlwr TYPE NO. I 

TEST CASE 

= ISflyO pound 
0 ? / i*?. PPOpELLANT 
OPTIMIZATION RUN 

SPECIFIC COv-^USTmn CHAVACTEhISTICS 


PER CFNT MtSS FUEL VAPOR I /FD 

FII 


100. 

PER CENT “ASS OF OxIOIZEM VAPORIZED 

F12 

C: 

100.000 

C* EFFICIENCY nOOEL 

ri3 


100.000 

C* EFFICIENCY PULSED CO^HtiSTORS 

FIA 


67.726A 

C* EFFICIENCY NO'l-PULSEO CONniJSTORS 

F15 



FUEL SYSTEM CHUOr-IN'G DECAY PATE 

F20 


-I26U. 

OZiniZFP SYSTE‘‘ CHUGGING DECAY PaTE 

F2I 


-I8I1. 

ptILSEO instability ChA° aCTERI ST I C 

FP2 

- 

-71.41G3 

NON-PULSEO instability CHftRACTEPlSTIC 

F23 


Not Cozcputcd 

pYKEKA FUEL STAPILIIY DFCAY RATE 

FZA 

- 

Not Computed 

PYKEHA OXIOIZEO stability OrCAY RATE 

F25 

£ 

-2920.49 

stability longitudinal TlvE LAG 

F2S 


-PP61.51 

stability TOAnSVEPSE TlfE LAG 

FP7 

C 

-P40.HP9 

STABILITY LRC PESPCNSE fUNCTlON 

F^A 

K 

-423.00ft 

stability Pf>IE“ LINEAR analysis 

F29 


-27334.7 

FUEL Pi^FSSuPE DPOP CHARACTERISTIC 

F31 


71 .25ftft 

pxtnr/hp p'>FbS’j*'r: r.-op c-'-^.'CTe- isT ic 

F32 

* 

SR-QflSA 

FUEL PLUS OXTD17ER HOLES Cmak ACTEM I STIC 

F41 

= 

S37.S14 

PXIOIZER dome volume CHARACTERISTIC 

FA2 

= 

10.3S73 

fuel dome volume C^*APACTEHlsriC 

F43 

= 

I0.5AB8 

PXIORFP hole LENGTH CmaPACTEm i ST I C 

F44 

«= 

•37ftS10 

FUEL mole length CmabaCTE^ISTIC 

fas 

- 

6.ftS3640E-02 

injector type Complexity characteristic 

FA6 

C 

Not Computed 

INJECTOR LF>*GTH ChAhactFRIST IC 

FSl 

c 

1 .9203% 

CHAmhFR. OlAHETrp CHArfACIt'wiSTIC 

F61 

= 

1.4B6S2 

MIXTURE RATIO CHARACTERISTIC 

F7I 


4. 94744 


TABLE 4-13. Specific Combustor Characteristics 
after 100 Design Perturbations 
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TABLE 4-14. Stability Root Set after 100 Perturbations 



Both roots are well damped; however, since the imaginary parts of 
these roots differ markedly from the value of Svh(3.0543) a "ray search" 
was carried out through the design space. This search proceeded along the 
ray joining the nominal engine design to the final design obtained after 
100 iterations. The ability to carry out this type of ray search through 
an n-dimensional space (in this case, a twelve-dimensional space) is a 
standard feature of the AESOP progreun. Fifty-two points were equi-distributed 
along the ray search joining the nominal and final design. The root corres- 
ponding to S-y}^ = 3.0543 was tracked along the ray starting from the nominal 
design. Root variation along the ray is presented in Figure 4-5. The root 
at 

z = (-.03257 + j3.0521) • ' 

presented in Table 4-7 tracks continually into the root at 

= (-.443094 + j2. 73775) 

confirming this root as a valid solution to the stability root characteristic 
equation. Both final roots, and Z2, obtained for = 3.0543 are, 
therefore, considered to be valid roots. Their heavily damped nature 
results in their providing no contribution to the final engine rating. 

It can be seen from Figure 4-5 that the root at Zj is becoming more stable 
as the design progresses and that the root at Z2 is becoming less stable. 

4.3.4 Verification of Optimal Solution . The optimal solution reported 
in Section 4.3.2 was verified in two ways. First jthe solution was continued 
for 100 additional iterations with the uniform random ray and pattern 
search algorithms. A slight performance improvement resulted. A final 
rating of 209.43 pounds was attained, a gain of 1.3 pounds over the solution 
of Section 4.3.2. Second^ the solution was restarted from the nominal 
solution using a different search algorithm. The algorithm used in this 
second solution was a recently developed directed rcindom ray search. 

Appendix B, in combination with the pattern acceleration algorithm. The 
final rating obtained by this method was 209.46 pounds after 100 iterations. 
Convergence of this solution has been added to Figure 4-1. It is clear 
that a final solution has been obtained. It is also clear that the newly 
developed search provides more rapid convergence to the solution than the 
older uniform directed ray search. This behavior is, in keeping with other 
tests of the new search. 


4.4 Conclusions for Sample Engine 

The AUTOCOM code has successfully developed an improved engine design 
starting from the existing nominal engine. The payload potential of the 
engine was improved by 87 pounds as measured by the rating equation. The 
computer time required by the AUTOCOM code was minimal. The average time 
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(8)1+ •* A oNmtauj 







requirement for an assessment of each combustor design wets approximately 
two seconds on the CDC 6600 computer. Computer time absorbed by the opti- 
mization subprogram AESOP in determining suitable design variable pertur- 
bations was negligible — approximately 103 seconds. The engine was optimized 
in one hundred design perturbations; hence, total computer time required to 
optimize the design was approximately four (4) minutes. More computer time 
would be required if combustor stability problems had been encountered. In 
this eventuality, it is estimated that twenty (20) minutes computer time 
would be required to obtain a solution. A definitive assessment of computer 
time in such a case awaits further experience using the AUTOCOM code. 

An examination of the optimal engine conponents reveals that the 
payload gain was largely obtained from improvements in the performance, 
chamber length, and chafer diameter characteristics. Small payload gains 
also resulted from improved stability and pressure drop characteristics. 

The injector complexity characteristic and the chamber mixture ratio charac- 
teristic both contributed performance losses when the final engine is 
conpared to the nominal engine. 

A complicated set of design variable perturbations were introduced to obtain 
the payload capability improvement. An assessment of the design variable 
changes by the optimization algorithms indicates that the number of fuel 
and oxidizer holes, volume of the oxidizer dome, volume of the fuel dome, 
length of the combustion chamber, chamber diameter, and the mixture ratio 
are all sensitive design variables in the engine considered. In particular, 
in both optimal solutions obtained the number of fuel and oxidizer holes 
rapidly rises to the upper limit permitted, indicating that further payload 
in 5 >rovement might result from a further increase in the number of holes 
allowed. Diameter of the fuel holes, diameter of the oxidizer holes, length 
of the fuel holes, and length of the oxidizer holes were relatively insen- 
sitive design variables for the engine design considered, presumably because 
of the basic stability of this engine. 
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B.l. INTRODUCTION TO THE OPTIMIZATION PROGRAM AESOP 


Program A£SOP (Automated Engineering and Scientific Optimization 
Program) is a digital computer subprogram designed for solution of a 
wide range of multivariable pai'ameter optimization problems on the 
CDC 6600, UNIVAC 1108, IBM 7094, or IBM 360 computers. The basic program 
described within this document has the ability to solve constrained non- 
linear optimization problems involving up to one hundred free parameters 
and twenty nonlinear functional equality or inequality constraints*. Of 
particular significance are the program's ability to 

1. Use a combination of search procedures ' 

2. Locate saddle points as well as maxima or minima 

3. Ix>cate more than one extremal 

Thirteen search techniques available for problem solution are listed in 
Table I. 'Phe searches may be employed separately or in any sequential 
combination, at the user's option. 

The AESOP program may be rapidly coupled to a wide class of parameter 
optimization problems including systems vdiich have been 5 >Tithesized as 
digital computer programs prior to their coupling to AESOP, For many 
problems this is achieved by adding a single named COMMON to the existing 
synthesis program and equating its controlling parameters to members of an 
array (ALPHAt^) find the output functions of interest to a second array 
(FUNCTNj). Large scale synthesis programs will require additional effort 
if the combined AESOP/synthesis programs exceed core size and hence require 
the use of the overlay feature. Treatment of this size of problem is 
enhanced in AESOP by the provision of a single call to the synthesis 
program from the optimizer main program. This call is then utilized by 
all searches permitting overlay of the. entire optimizer against the 
synthesis program. 

The optimization program is v/ritten in the algebraic . FORTRAN IV 
language. Tne program has been successfully operated on the CDC 6600, the 
UinVAC 1108, the IBM 7094, and the IBM 360. 

The subprogram AESOP may readily be retrieved from the AUTOCOM program 
and treated as a subprogram entity. 
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TABLE I. DESCRIPTION OF THIRTEEN SEARCH 
METHODS 
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MAPPING SBARQi : Generates a contour map of the performance surface perturbing two parameters at a time 

holding all others constant. 


B.2. A GUIDE TO PROGRAM AESOP USAGE 


B.2.1 Outline of Program AESOP 


AESOP is constructed in a modular fashion. It includes more than 
one hundred subroutines which, for descriptive purposes, can be arranged 
into the thirteen sets illustrated schematically in Figure 1. Each set 
of subroutines collectively function to perform one of the searches 
described in Table I. In addition to the thirteen major subroutine sets 
various ancillary subroutines are used in support of the search procedures. 
These include, for example, the multiple extremal search procedure and the 
saddle point search procedure. 

. A flexible input format is provided through the FORTRAN N,y>ELIST 
input procedure. All AESOP data is entered through a single NANELIST 
block labelled lAESOP. All optimization data is nominally set by block 
data within AESOP; this data may be modified at the user's option; only 
that data which varies from the standard value need bev read into the 
program. This feature reduces data input and hence the tedium, of preparing 
a case or the chance of input error. 

The program structure employed in the synthesis program is entirely 
dependent on the particular problem being studied. For example, in 
Reference l,an applicatipn to the design of hypersonic cruise vehicles, 
the MAIN program controlled an extensive overlayed program containing all 
major technologies entering into the design of a high speed aircraft. The 
structure for the polynomial valley function (j) = fax - b)^ + (cy - d)N = exy,. 
on the other hand, would consist of a single coded expression.. The opti- 
mizer structure is independent of the synthesis model employed.' The 
optimizer selectively calls the thirteen search controlling subroutines. 

Each one of the thirteen search controlling subroutines may split into 
several lower level subroutines controlling major searcli functions. The 
two auxiliary subprograms WARPS (for multiple extremal' search) and SADDLE 
(for solution of mini-max problems) may . be used in , conjunction, with any 
search procedure or search procedure combination. 

More detailed information is presented in Reference 2. Detailed 
results illustrating the application of program AES.OP are contained in 
References 3,4 ,5 ,6 , and 7. The analytical development of the 
program is contained in References 2 and 3. 


B.2.2 Coupling of Mathematical Model to the Optimizer 

Program AESOP assumes that it is operating on a programmed model of 
a nonlinear parametrically defined system and that a system evaluation 
with a given set of parameter values results in a unique evaluation of 
the system output functions. 
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Figure 1. Schematic of Optimization Program - AESOP 


















Several modifications are required in the system main program, 
specifically: 

1. COMMON/ AE SOP D/ , of the optimization program 

should be introduced. 

2. The parameters to be optimized should be equated or 
equivalenced (being careful to avoid difficulties 
with COM^^ON) to the ALPHA array contained in COMMON/ 

AESOPD/. 

3. The functions of interest' which include 'the perfor- 
mance and constraint functions should be equated or 
equivalenced to the array FUNCTN contained in 
COMMON/ AESOPD/ . 

NOTE: The performance function must be phrased as 

a function to be minimized. 

4. Any system input data in addition to the optimization 

program input should only be read when JJJ = 1 . This 
integer variable contains the number of performance 
evaluations performed by the optimizer and is also 
contained in CONE-ION/AESOPD/ . ’ ■; 

5. An initial call to AESOP must precede the optimization 
loop. The optimization loop must be created using 
the integer NSR contained in 'COMMON/ AESOPD/. 

This can be illustrated ^^?ith a simple example. S'uppbse the function 
<|) = lOO(y-x^)^' + (1 - x)2 is to be minimized. A satisfactory CDC 6600 
MAIN program would be 


-PROGRAn MAIN -CINPUTv OUTPUT-, TAPES=INPUT-, TAPEb=0UTPUT-, 
$ TAPE1-. TAPES} 

COnnON /AESOPD/ ADATA-CSODD} 

C DEFINE AESOP COriNON 

REAL ALPHA-QOD} 

C DEFINE THE PARAI1ETER ARRAY SIZE 

REAL FUNCTNaaO} 

C DEFINE THE SYSTEN FUNCTION ARRAY SIZE 
ECUIVALENCECADATA-C 7^a}. ALPHA } 

C DEFINE POSITION OF FIRST PARAMETER 

ECUIVALENCE-L ADATA-CSSSS} ^ FUNCTN} 
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C DEFINE POSITION OF FIRST FUNCTION 

' E(3UIVALENCE-CADATAOL2S}-. NSR > 

C DEFINE POSITION OF AESOP LOOP INDICATOR 

EtJUIVALENCEIALPHAClI-. XlnCALPHACaii YI 

C E(3UATE SYSTEH PARAHETERS TO AESOP OPTiniZATION PARAMETERS 
CALL AESOP 

C READ IN ANY NON STANDARD OPTIMIZATION DATA 

7777 CONTINUE 

C PARAMETER OPTIMIZATION LOOP 

F = 1DD.0*{Y -X**2>**E +-C1.0 -X}**2 

C EVALUATE SYSTEM 

FUNCTN-C1}=F 

C TRANSFER SYStEM OUTPUT INTO AESOP FUNCTION 

CALL AESOP 

C AESOP OPTIMIZATION POINT 

IF {NSR .E(3. 1} GO TO 7777 

C RETURN FOR NEU EVALUATION UNTIL SYSTEM OPTIMIZED 

END 


B.2.3 AESOP Overlay Structure 

For programs that., exceed core size the overlay structure shown in 
Figure 2 may be used. For the overlay version of AESOP it will be 
necessary to make several changes to the MAIN PROGRAM for links 0,0 and - 
1,0; the setup for overlay is shown in Figures 3 and 4. Figures 
presents the overlay version of SELECT, and Figure 6 shows the additional 
programs required by the overlay structure. The text cards for COPYLIB 
are given in Figure 7 The total core required to load the AESOP over- 
lay configui'ation showi in Figure 2 is approximately SOOOOg. 


97 






oe: 



UJ 

U3 

u 

Ui 

CO 


c: 

29 


►H 

< 

u 


*-H 

X 

tj 

t; 

Z3 


< 

to 




o X J J 

a. Lu <n tu < o 
tuc!:coHu: 2 Hf-<ca 
>>cno<;ujoos: 
Ocitu.JcxE-^>-l>- 
SQSOuE-XCuCJucn 


/-><iJ 5 ci 03 — :<H 
o><u^sc/:cj 2:2 


' H >- 


CM § S ^ S b tD G u CU 

^u.oocooocucL.r: 


r-i 


a 

Z S j tx z ^ u 

0 

Ctf H-* < H -3 cc 2: >- 


HJ>UJ<UOUU 

rC 

totxuco>E-«G;i-H^ 


w J z ui < < to tu z 


oo<i-Htxcoa cxstu 


e, fe o 


X ex cj >- 

to < o •• f- 2 : 

coJf-ieicoi— « 02 ;<s; 

tl 1 tu: 2 <UJ<<U 4 U 

<;to>-» 3 ccQs:tocx£u 


Kl C£ Dl. G- 

e: ir ^ 23 

M« A* 

•U ex H 2; 
r>j ei < !d O 
'^to A O U 


(2,14) 

SRCH12 

JACOBS 

GRADV 

(2,13) 

SRCHll 

RAYSEC 

^0 >- CO 










fvl r -4 < 0 

> 









— 1 X K u 

< z 









^0 z 5 ^ 

CO UJ 









CM C£ < 

m 0 









<-^to a c 2 

cx cs 









/-Ncn >- to 










0 < 0 

>■ 









r-l X C 2 CJ 

< z 









•*CJ Z 2: 

to u 









fM ei < < 

CU CJ 









^(0 K CC 

CX CS 









#^oo H 










0 0 z 

> 









X ►-* Z 

< 









. 0 0 UJ 

CO 









rsi ci ex 0 

UJ 









^ CO cc oi 

CX 









z 

X Q z 

CsS 

X 







0 Q > 

es H C w 

X 

>- 






:rf HH ►— ( 

H CO H 

X 








•vcj > os 

< t-H CJ 

> 

d 

z 






CN) os < Ul 

^ z u 

< 

X 

X 






^to Q Q 

Q hH to 

to 

CO 

X 






/-^ \o <C ex 

H -J Z 

a. 

X 

X 

z 


> 

X 

X 

vo 0 ex — 3 

CS X < 

0 X 

X 0 


z 

X 

>- 

- X a < 

UJ Z OS 



to 



»— ♦ 

< 

X 

CM CJ < z 

ex ix f— • 

< 

X 


b 


X 

X 

z 

w CS X t-H 

XX a 

X 

< z 

X 


Zil 


X 

CO O S 

< < c/cy 



CO 

z 

O' CO 

X 

to oS X 

x; 









10 0 ex X 

X > 

X 








X tU Q 

CJ < z 

X 








<N tj tU Z 

> to X 

X 








cS OS <L 

CO X 0 

“ 








to U K 

< U, CS 

CO 








CJ 

Q Z 

CX 

Cu 







^ 0 CO > 

H 0 

X 

>- 







^ X LJ hh 

H to H 

< 

X 







<sj CJ Q OS 

< *-t CJ 

X 

z 

i 






w eS H u 

^ z X 

H 

X 







CO CO Q 

Z CO 

to 

X 

1 



. 



to >* 










^ 0 fx 










to X 2; 










^ t_J 0 










rvi OS <; 










^ CO 










/ — \ G. 00 

< CO X 

CJ 

£0 

X 

X 

X 




CM H 0 < 

H H t-- 

f—, 


X 

X 

X 




-* X CO E- 

X X tD 

S— . 

X 

X 


X 




CM ex 01 <: 

000 

0 

0 

0 

0 

0 




Z < Q 

XXX 

X 


X 


X 




i-i 03 CO 

z z z 

z 

z 

z 

z 

z 





ST g« tx 

i-HoorDiuo -JxJ 

cMtju >;:duixg: 2 cq 
wctcuj-^Og:o:::: 3 ^ 

LOtOCSCOtOCitOtXtO 


e> 

CO o 
•ri 
CO 

^ ti) -w 

o ^ ;s 
•^ SJ o 
^ s: X 

10 3 

to 


98 


Figure 2. AESOP Overlay Structure for CDC 6600 







PROGRAn nAIN{INPUT-, OUTPUT} n TAPES=INPUTt TAPEb=OUTPUTi 
, $ TAPE1-. TAPE5} 

connoN/AESop])/ ADATATanao} 

El3UIVALENCE-CAI>ATA-C31E5}-. NSR 
100 CONTINUE 

CALL OVERLAY-CTLAESOPOLi 2^ Oi LHRECALLl 
7777 CONTINUE 

CALL OVERLAYCTL'AESOPOL-. It 0-. LHRECALLI 
CALL 0VERLAYC7LAES0P0L-> 2-. 0-. LHRECALL} 

IF -CNSR -EC2. 1} GO TO 7777 
GO TO 100. _ 

END ' ' 


Figure 3- Sample MAIN Program for Link 0, 0 


PROGRAM Et2NS 

COmON/ksOPO/ AOATA-CaOOD} 

• REAL ALPHA-CIOOI 
REAL FUNCTNODDI 

EfiUIVALENCEIADATAI 7M2>. ALPHA }. 
Et3UIVALENCEIAI)ATAi2SS2}-i FUNCTN} 
EQUIVALENCE! ALPHA-Clln XItCALPHA-CE}-. YI 
C ROZENSROCK VALLY . 

FUNCTN-CII = 100.D:«:-CY -X2>*:«:2 +-C1.0 -X}**2 
■ END 


Figure 4. Sample MAIN Program for Synthesis Link 
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SUBROUTINE SELECT 
COMMON /AESOPD/ADATA(8000;) 

EQUIVALENCE (ADATA(2952) , METHOD) 

GO TO (1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11, 12, 13), METHOD 

1 CONTINUE 

CALL OVERLAY (7LAESOPOL, 2, 1, 6HRECALL) 

RETURN 

2 CONTINUE 
RETURN 

3 CONTINUE 

CALL OVERLAY (7LAESOPOL, 2, 3, 6HRECALL) 

RETURN 

4 CONTINUE 

CALL OVERLAY C7LAESOPOL, 2,4, 6HRECALL) 

RETURN 

5 CONTINUE 

CALL OVERLAY (7LAESOPOL, 2,5, 6HRECALL) 

RETURN 

6 CONTINUE 

CALL OVERLAY (7LAESOPOL, 2,6, 6HRECALL) 

RETURN 

7 CONTINUE 

CALL OVERLAY (7LAESOPOL, 2, 7, 6HRECALL) 

RETURN 

8 CONTINUE 

CALL OVERLAY ( 7 LAESOPOL, 2, lOB, 6HRECALL) 

RETURN 

9 CONTINUE 

CALL OVERLAY(7LAESOPOL, 2, IIB, 6HRECALL) 

RETURN 

10 CONTINUE 

CALL OVERLAY (7 LAESOPOL, 2, 12B, 6HRECALL) 

RETURN 

11 CONTINUE 

CALL OVERLAY C7LAESOPOL, 2, 13B, 6HRECALL) , 

RETURN 

12 CONTINUE 

CALL OVERLAY(7LAESOPOL, 2, 14B, 6HRECALL) 

RETURN 

13 CONTINUE 

CALL OVERLAY (7LAESOPOL, 2, 15B, 6HRECALL) 

RETURN 

ENTRY AESOP I 

CALL OVERLAY (7 LAESOPOL, 2,2, 6HRECALL) 

RETURN 

END 


Figure 5. Overlay Version of Subroutine SELECT 
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J 


Figure 


PROGR.-y-{ SRCH20 
CALL AESOP 
END 

PROGR,A^f SRCHOl 
CALL SECCON 
END 

PROGRAM INPUTA 
CALL B AESOP 
END 

PROGR.AiM SRCH03 . 
CALL MAGNFY 
END 

PROGRAM SRCH04 
CALL STDESC 
END 

PROGRAM SRCH05 
CALL CREEPR 
END 

PROGRAM SRCH06 
CALL QLFADRA 
END 

PROGRAM SRCH07 
CALL DAVIDN 
END 

PROGRAM SRCH08 
CALL RPOINT 
END 

PROGR.A^^ SRCHO'9 
CALL R.ANRAY 
END 

PROGRAM SRCHIO 
CALL RAXRAY 
END 

PROGR.AM SRCHll 
CALL RAYSEC 
END 

PROGRAM SRCH12 
CALL JACOBS 
END 

PROGRA! SRGH13 
CALL MAPPER 
END 


Additional Programs Required by the Overlay Version 
of SELECT 
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OVERLAY (AESOPOL, 0,0) 
MAIN 

OVERLAY (AESOPOL, 1,0) 
EQNS 

OVERLAY (AESOPOL, 2,0) 
SRCH20 

OVERLAY(AESOPOL,2,l) 

SRCHOl 

OVERLAY (AESOPOL, 2, 2) 
INPUTA 

OVERLAY (AESOPOL ,2,3) 
SRCH03 

OVERLAY (AESOPOL, 2, 4) 
SRCH04 

OVERLAY (AESOPOL, 2, 5) 
SRCH05 

OVERLAY(AESOPOL,2,6) 

SRCH06 

OVERLAY(AESOPOL,2,7) 

SRCH07 

OVERLAY(AESOPOL,2,10) 

SRCH08 

OVERLAY (AESOPOL, 2, 11) 
SRCH09 

OVERLAY (AESOPOL, 2, 12) 
SRCHIO 

OVERLAY (AESOPOL, 2 , 13) 
SRCKll 

OVERLAY (AESOPOL ,2,14) 
SRCH12 

OVERLAY (AESOPOL, 2, 15.) , 
SRCH13 


Figure 7. COPYLIB Text Cards 
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B.3. AESOP DATA INPUT 


Program AESOP data can be conveniently grouped according to search 
and function. The user employing a particular search can independently 
specify the characteristics of that search and may not be concerned with 
input relevant to the other searches. Hence, a data grouping by search 
and by function is presented below for user ' convenience . It should be 
noted that certain inputs are common to more than one search; where this 
occurs, the input is repetitively defined in each search. All AESOP data 
is introduced through a single FORTRAN NAMELIST labelled lAESOP. All 
data is nominally established by the Block Data Subroutine BDATA7. 


B.3.1 Search Selection and Control 

There are thirteen search procedures embodied in the AESOP subprogram. 
The user may select any sequential combination of these searches through 
data input. 


^^•CMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

NUMOPT 

The number of optimization techniques to be 
employed. Each individual search request in 
a sequence of requests adds to this input 
(e.g., the search sequence 4,2,d,2 requires 
NUM0PT=4). Maximum number of searches 
employed must satisfy NUMOPT < 20. 

2 

METHOP- 

The search sequence by. numeric identification. 

For example, the input METHOP(l) = 1,2, 3, 4, 5, 
6,7,8,9,10,11,12,15 signifies t]\e following 
search sequence: 

1 - Sectioning 

2 •* Pattern 

3 “ Magnification 

4 - Steepest-Descent 

5 - Adaptive Creeping 

6 - Quadratic 

7 • Davidon (Fletcher-Powel 1 ) 

8 “ Random Point (Montc-Cario) 

9 - Random Ray (random evolution) 

10 - Directed Random Ray 

11 “ Arbitrary Ray 

12 “ Jacobson j 

13 ' Mapper 

10, 2 


The complete search sequence will be referred i 

to as an ovtimizdtion cijcle. 


mxjjj 

1 

The maximum number. of system evaluations, A 
direct iteration ynufzber • limit . 

The maximum number of times the search sequence 
(optimization cycle) defined in MI:TH0Pi will 
be utilized. 

200 

I RE PET 

1 
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B.3.2 Parameter Selection 


The user may equate any parameter in the system under study to an 
optimization parameter. This is usually most readily accomplished in 
the user's MMN program. If the initial AESOP call follows the system 
program, input parameter values entered through AESOP override the 
system values. 


MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

NALPHA 

The number of parameters available for opti- 
mization. No more than one hundred param- 
eters may be employed. 

0 

ALPLOi 

Lower bounds on each parameter search range. 

100*0.0 

ALPHIi 

Upper bounds on each parameter search range. 

100*0.0 

ALPHAi 

The nominal parameter values. Note that 
ALPLO^^ALPHAifiALPHIj^ must be satisfied. If 
a particular parameter, say ALPHA-; , is to be 
fixed in value in a particular computation, 
then set ALPLOj=-ALPlLAj=ALPHIi . This effec- 
tively reduces the parameter space dimension 
by one for each such parameter, 

100*0.0 

HILO 

If HILO is greater than zero the lower and 
upper bounds on each parameter will be 
dynamically determined such that 

HILO* 1 ail 

OLi = oti 

. HILO* 1 ail 
«hi = * 2 

provided HILO * \<^i\ > (ct*hi ~ where 

ct^Li Lhe current values of the 

parameter bounds. This option permits the 
1 user to define dynamically a percentage 
1 search region centered about the current 
! value for each parameter. 

0 
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B.3.3 Multiple Extremal Option 


, AESOP contains a multiple extremal search procedure based on a 
warping of the response surface, _ In this transformation hypervolumes 
in the region of a knoim extremal .point are reduced in magnitude and 
hypervolumes far removed are increased in magnitude. 


MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

I WARP. 

Controls multiple extremal option; 

= 1, automatically warps the response 
surface 

= 0, leaves the response surface unmodified 

0 

WARPALi 

The point at which the warping transformation 
is centered, i.e., the location of a known 
extremal point 

100*0.0 

WARPN 

The degree of the warping transformation. The 
greater WARPN, the greater the response 
surface distortion. 

2.0 


B.3.4 Optimization Function Selection 


AESOP operates on an internally defined array of functions. The 
function to be minimized and those to be constrained are defined by 
data input. This permits ready specification of alternate optimization 
problems. It should be noted that constraints are imposed on the basis 
of an external penalty function approach using quadratic penalty function 
terms. Penalty function weights are adaptively determined within the 
code. 


MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

FUNCTNi 

NFUNC 

An internal array containing all computed 
optimization functions 

The total number of functions (FUNCTf'Ji) being 
defined and computed in the system model. 
NOTE: NFU^^C 100. 

100*0.0 

1 
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MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

NPHIAC 

The function to be minimized. AESOP always 
searches for a minimum; to maximize FUNCTNjn 
define 

FUNCTN = -FUNCTN' 

, . . . n ^ 

and minimize 

FUNCTNn 

1 

NUMPSI 

The total number of functions being con- 
strained. Note that NUMPSI <20. 

0 

NPSI^ 

The functions to be constrained, e.g., NPSI(l) 
equal to 3 , 5 , 1 , 7 indicates that FUNCTN3, 
FUNCTNv , FUNCTNi, and FUNCTN7 con- 

strained. 

20*0 

SIBARi 

The desired values of the constraint functions 
defined by NPSIj^. 

20*0.0 

FTOLi 

The acceptable final tolerances on the con- 
straint function values, SIBARj^. 

20*1.0 

TTOLi 

Initial acceptable tolerances on the constraint 
function values (siiould be approximately 100 
times greater than the corresponding FTOLi)-. 
These tolerance bands are subsequently reduced 
in size at a problem dependent rate until they 
equal the appropriate FTOLi* 

20*100.0 

PSIlVTi 

Initial constraint error weighting factors in 
the augmented performance function, (()*, where; 

i 

20*li0 


= + 1 WiCt;;. - ii;.)2' 

-i- ■ . , , 



Here = PSIlVTi 


WTUPi 

Incremental multiplicative constants used to 
increase’ the oh constraints which prove 

difficult to satisfy. The' nominal values of 
WTUPi= 2.0 should be acceptable. 

20*2.0 


< 
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MNEMONIC 

DESCRIPTICN 

— 

NOMINAL 

VALUES 

WTDOi^'i 

Decremental;inultiplicative constants used to 
decrease the when a constraint is easily 

satisfied. The nominal values of WTDO\VNi=0.5 
should be acceptable; hence, this input can 
normally be omitted. 

20*0.5 

TOLFAGj^ 

, Constraint tolerance reduction factor. 

20*0.5 

INDPSI 

Constraint indicator; if INDPSI 
>0, will be constrained 

<0,- will be constrained <\\)^ 

«0, will be constrained to be within 

±FT0L of ’J;. . 

Yl 

0 

IPDWGT 

Print indicator for constraint weights. 
Equals 1 -to obtain print out of constraint 
weights each time they are modified. 

0 


B.3.5 Saddle Point Search 


Saddle points are located by replacing a function minimization 
criteria. At the user' s . option a sign of curvature correction may be 
introduced as discussed in Reference 5, 


MNEMONIC 

. DESCRIPTION 

NOMINAL 

VALUES 


Used to indicate minimizing or maximizing 
control for each parameter 
-1 indicates minimizing control 
+1 indicates maximizing control 

100*-1 

SPOWER 

Exponential power used in computing the 
penalty term when the curvature is the 
wrong way on a given control parameter. 

10.0 

SFACTR 

Control step size for computing the gra- 
dient (i.e., .05 will produce a perturba- 
tion of 5% of the range 

1 

0.05 
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MNEMONIC 

DESCRIPTION 

NOMINAL 
VALUES . 

ISADEL 

=0, saddle point search not used 
=1, saddle solution only 
=2, minimuTn only 
=3, maximum only 

=4, saddle, minimum or maximum (minimum 
gradient magni tude) 

0. 

I PAS AD 

Print indicator for alpha array; equal to 
1 to obtain alpha print from saddle 

0 : 

IPFSAD 

Print indicator for function array; equal 
to 1 to obtain function print from saddle 

O' 

IPDSAD 

Detail print indicator; equal to 1 to 
obtain detail print from saddle 

0 


B.3.6 Detailed Search Data 

Data for each of the ten searches is defined below. 

B.3.6.1 S ectioning Search Data (METHOPj = 1) . Uses the Golden Section 
method to carry out a sequence of one parameter at a time optimization 
searches. 


' 

MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

LIMIT 

Tne number of times each parameter will be 
sectioned during a single request for 
sectioning search . 

2 

NMAXLO 

Maximum number of point evaluations employed 
in a single parameter's sectioning at search 
commencement (first optimization cycle) . In 
successive cycles, the maximum number of 
points employed is increased by one. 

15 , 

ISECOF 

Indicator used to '^tum off” the sectioning 

1000 

1 

search. \Vhen the number of optimization 
cycles exceeds ISECOF, the sectioning 
search will be ”tumed off.” 
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MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

PHIEPS 

Performance values within PHIEPS of the 
minimum value yet attained are treated as 
being equal in the Golden Section 

0.0 

NMAXUP 

An upper bound on the maximum number of 
point evaluations employed in sectioning 
a particular parameter 

25 

ISIDE 

Indicator specifying an arbitrary selection 
of left or right boundary for a parameter 
that does not appear to affect the system 
performance 

ISIDE=0, select lower limits 
ISIDE=1, select upper limits 

0 

XTENHI. 

1 

Optional extension of higher search limits 
(ALPHIj^) for a parameter that has gone to 
the ALPHI^ boundary 

100*0,0 

XTENLOi 

Optional extension of lower search limits 
(ALPLOj^) for a parameter that has gone to 
the ALPLO^ boundary 

100*0.0 

IRANDM 

. . . . 

j 

1 

Controls the order in which the parameters 
are sectioned: 

= 0, random order selected 

=1, natural order selected 

=2, reverse natural order selected 

=3, order of selection is input via KORDER 

0 

KORDER 

An array that contains the subscripts of the 
parameters in the order that they are to be 
selected. Required input if IRANDM=3; other- 
wise generated automatically by AESOP 

100*0.0 

FACTHI 

Section termination criteria. If three 

f 

successive performance function values are 
within FACTHI of each other during 
sectioning of a given parameter on the 
first optimization cycle.,, the section search 
i of that parameter will cease. The termina- 
tion criteria is internally reduced with 
each optimization cycle. 

l.OE-3 
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MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

FACTLO 

The lower limit, on the- termination criteria 
in any optimization cycle 

l.OE-5 

ITRADE 

Optimization/trade study indicator: 

= 0 , carry out a normal optimization search 
= 1 ,,^ determine performance function .sensi- 
tivity to each parameter by sectioning 
each parameter in turn about a given 
fixed point in parameter space. 

0 

IPASEC 

Print indicators for parameter, detail. 

a 

IPDSEC 

and function values during sectioning 

0 

IPFSEC 

search. : . 

.0 


B.3.6. 2 Pattern Search Data (METHOP-j = 2) . ‘Pattern search is a 
simple generalized acceleration procedure which may be employed to 
attempt an acceleration of the trend revealed by a previous search 


MNEMONIC 

> DESCRIPTION 

, NOMINAL 
VALUES 

I PAP AT 

Pattern search print indicators for param- 

0 

IP DP AT 

eter detail and function values during 

0 

IPFPAT 

pattern search. 

0 

NMAXUP 

Maximum number of pattern steps;, otherwise 
search will end when cjinew > <f>old; i.e., 
failure to "gain" performance 

.25 

B. 3.6.3 

Ma 0 iification Search Data (METHOP^ = 31. 

This search 

c . 


useful when the unconstrained extremal solution is the null vector. It 
searches the ray joining the origin to the current point in parameter 
space. , . 


MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

MAXMAG 

Maximum number of point evaluations per- 
formed during a single magnification 
search ’ * 

99 
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MNENOMIC 

DESCRIPTION 

NOMINAL 

- 


VALUES 

DELMAG 

The magnification perturbation size, nom- 
inally set to one per cent of distance to 
origin. Not normally modified from nominal 
value. 

0.01 

IPAMAG 

Point indicators for parameter, detail, and 

0 

IPDMAG 

function values during each magnification 

0 

IPFMAG 

k-a 

search . ' 

0 


: B.3.6.4 Steepest-Descent Search Data (METHOP^ ' = 4) . The AESOP 
steepest'descent search features several metric tensor options to 
deflect the search in a more promising direction than" the gradient 


MNEMONIC . 

.DESCRIPTION 1 ... 


NUMSTD 

Number of gradient evaluations and one- ' 
•dimensional searches performed each time 
that a steepest-descent search is 
requested during the optimization cycle 

2 

i'ndwma 

Steepest-descent weighting matrix indi- 
cator 

=0. , unit matrix 
= 1 . , eriipirical matrix ' / 

=2., alternate on each cycle between unit 

0 

■ 1 


and empirical* matrices 


WITER^ 

Learning factors for steepest-descent 
weighting matrix 

100*1.0 

NMAXLO ' 

Maximum number' of point evaluations employed 
in the steepest-descent one-dimensional 
ray search at search commencement (first . 
optimization cycle). In successive cycles 
the maximum number of point evaluations 

15 


permitted is . increased by ,1 . .. . , . . 


NMAXUP 

Upper bound on the number of point evalu- 
ations along a steepest-descent one-dimen- 
sional ray in any optimization cycle. 

25 
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MNEMONIC 

DESCRIPTION 

NOMINAL 
VALUES • 

FACTHI 

One- dimensional steepest-descent ray search 
termination criteria during first cycle. The 
termination criteria is reduced in each 
successive optimization cycle 

-l.OE-3 - 

FACTLO 

Lower limit on one -dimensional steepest- 
descent ray search termination criteria, -in 
any optimization cycle. 

l.OE-5 

XTENHIi 

Optional extension of higher search limits 
(ALPHI^) for a parameter that has gone to 
the upper bound 

100*0.0 

XTENLOi 

Optional extension of lower search limits 
(ALPLOj^) for a parameter that has gone to ; 
the lower bound . 

100*0.0 

PHIEPS 

Performance values within PHIEPS pf the 
minimum value yet attained are treated as 

0.0- k; 


being equal in the Golden Section 

, 

, DALFMN 

Minimum parameter perturbation ■ . • ; 

5.0E-4 

IPW^T , 

Weighting matrix print control indicatory 
= 1, to print weighting matrix ■ ; 

0 

IPDRIV 

Controls print out of derivatives:. ; • 

=1, to print derivatives 

0 

IPASTD 

Steepest-descent search print, control indi- 

' 0 . . 

IPDSTD 

cators for parameter, detail, and function 

0 

IPFSTD 

values. - 

0 


B.3.6.5 Adaptive Creeping Search Data (MHTHOPi = 5) . An elementary 
univariate search procedure which arrives at an optimum by a. sequence 
of small parameter value perturbations. Perturbation size and direction 
are adaptively determined within the code. 
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DESCRIPTION 

' 

NOMINAL 

VALUES 

MAXCRP 

Number of creeping search perturbations 
introduced into each parameter by a single 
adaptive creeping search in the optimization 
cycle. 

5 

IRANDM 

Controls' the order in which parameters are 

perturbed: 

rO, random order 

=1, natural order 

=2, reverse natural order 

=3, order of selection is input via KORDER 

0 ■ 

KORDERi 

An array that contains the subscripts of 
the parameters in the order- that they are to 
be selected. Required input if IRANDM=3. 

100*0 

DCREEP^ 

The initial perturbations to each parameter 

100*1. OE-3 

CREPMNi 

Minimum perturbations for each parameter 

100*1. OE-7 

CREPMX^ 

Maximum perturbations for each parameter 

100*1. OE+6 

ALFSINi 

Direction of perturbation for each param- 
eter (ALFSINj=±1.0) . May be used when 
picking up a search sequence to indicate 
the improving direction 

ioo*i!o ■ 

IPACRP 

Adaptive creeping search parameter, detail 

0 

IPDCRP 

and function values print indicator 

0 

IPFCRP 


0. , 


B.3.6.6 Quadratic Search Data (METHOPt =6). The quadratic search 
is a second-order search procedure formally identical to the Newton- 
Raphson search. 
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MNEMONIC 

DESCRIPTION 

NOMINAL • 
VALUES 

QPERTi 

Parameter perturbation magnitudes employed 
in computation of numerical partial deriv- 
ative matrices 

1 . . , 

9^4) , d(j> ' 

dajdttj daj^ 

20*5. OE- 3 

QFACTR 

Scaling factor on the QPERT^^ 

1.0- - 

NMAXLO 

Maximum number of point evaluations employed 
in the quadratic one-dirnensional ray search 
commencement (first optimization cycle) . In 
successive optimization cycles, the number 
of point evaluations permitted increases by 
one. 

15 

NNL4XUP 

Upper bound on the number of point evalua- 
tions along a quadratic one-dimensional ray 
search, in any cycle. 

25 

FACTHI 

One-dimensional quadratic ray search . termi- 
nation criteria during first cycle. The 
termination criteria is decreased in each 
i successive optimization cycle. 

l.OE-3 . 

FACTLO 

Lower limit on one-dimensional quadratic 
ray search termination criteria, ir\.any . . 

cycle 

l.OE-5 

IP AQUA 

Quadratic search parameter, detail, and 

0 

IPDQUA 

function values print indicators 

■ 

0 • 

IPFQUA 

0 


B.3.6.7 Davidon Search Data (^■lETHOP^ = 7) . Davidon's method is a 
pseudo-second order steepest-descent method. Second order information 
regarding the response surface characteristics is obtained indirectly 
and incorporated in the metric tensor. 


114 






MNEMONIC,. 

DESCRIPTION 

NOMINAL 

VALUES 

MAXDVD 

1 

.Number of Davidon (Fletcher-Powell) gra- 
dient evaluations and one-dimensional 
searches performed each time that a Davidon 
search is requested in the optimization 
cycle 

10 

QPERTi 

Parameter perturbation magnitudes employed 
in computation in numerical partial deriv- 
atives , 

; 

100*5. OE- 3 

NMAXLO 

Maximum number of point evaluations employed 
in the Davidon one-dimensional ray search at 
commencement (first optimization cycle). In 
successive optimization cycles, the number 
of point evaluations permitted is increased 
by one 

15 

NMAXUP 

Upper bound on the number of point evalua^ 
tions along a Davidon search one-diinensional 
array in any cycle ' . , 

25 

FACTHI 

One-dimensional Davidon ray search termina- 
tion criteria during first optimization 
cycle. The termination criteria is 
decreased on each successive optimization 
cycle 

! 

l.OE-3 

FACTLO 

bower limit on one -dimensional Davidon 
ray search termination criteria, in any 
optimization cycle 

l.OE-5 

DALFMN 

IPDRIV 

Minimum parameter perturbation 

Controls print out of derivatives: 
to print derivatives 

5.QE-4 

IPADVD 

Davidon search parameter, detail and 

O ' 

IPDDVD 

function values print control indicators 

0 

IPFDVD 

• 

0 
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B. 3.6.8 Random Point Search (NETHOP^ = 8) . A straightforward 
Monte-Carlo procedure involving a uniform distribution of random 
points in the feasible region. The search automatically turns ‘ 
itself *loff*^ after the first optimization cycle. 


MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

MAXRPT 

The maximum number of random points to be 
employed in the first request for a random 
point search within the optimization cycle. 
In successive requests, MAXRPT is set to 
zero and no evaluations result. 

10 

IPARPT 

Random point parameter, detail, and function 

0 

IPDRPT 

values print indicators. 

0 

IPFRPT 


0 


B.3.6.9 Random Ray Search (METHOPj^ - 9) . An optimal point is found 
by means of an evolutionary sequence of small small random steps 
in the control space. If a random step is unsuccessful, the antiparallel 
step is always introduced. 


MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

MAXRRS 

The maximum number of random rays , one or 
two sided, investigated each time the opti- 
mization cycle requests a random ray 
search 

3 

RUFHI 

The initial maximum non-dimension pertur- 
bation measure for each parameter. This 
is reduced each time random ray search 
consistently fails to improve performance ^ 

I.OE-2 

RUFLO 

Minimum-maximum dimensional perturbation 
measure for each parameter 

I.0E-I3 
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MNEMONIC 


DESCRIPTION 


NOMINAL 

VALUES' 


SPRANX Starting trigger for the random number 33333333333 

generator 

RUDOIVN ' Used to adaptively' alter (increase or .25 

decrease) the step size for the random 
ray search 

INGRAN Number of sequential two sided unsuccess- 6 

ful random steps required before the step 
size will be reduced 

(AS = AS * RUDOltfN) 

lUPRAN Number of sequential two sided successful 6 

random steps required before the step size 
will be increased 

(AS = (AS/RUDOM + AS)/2.0) 

RANGEN Indicator for random number generator 1.0 

>0, generates imiformly distributed random 
numbers 

<0, generates normally distributed random 
numbers 

FNORMAL Value by which the normally distributed 5.0 

random numbers will be divided in order to 
reduce the 3a value of the range over which 
random numbers will be generated. If 
FNORMAL = 1, the 5a value will be approxi- 
mately ±Tr 

IPARRS Random ray search parameter, detail, and 

IPDRRS function values print indicators. 

IPFRRS 


B. 3.6. 10 Directed Random Ray Search (METHOPj = 10) . This search 
proceeds along a succession of random rays distributed about a 
best estimate of the gradient vector. 


117 


o o o 






I^OMINAL 

VALUES 


MNEMONIC 

DESCRIPTION 

I^OMINAL 

VALUES 

MAXRRS 

The maximum number of random rays, one or 
two sided , investigated each time the opti- 
mization cycle requests a random ray 
search 

3 

. 

RNFHI 

The initial maximum non-dimension pertur- 
bation measure for each parameter. This 
is reduced each time random ray search 
consistently fails to improve performance. 

l.OE-2 

RNFID 

Minimum-maximum dimensional perturbation 
measure for each parameter 

l.OE-13, 

SPRANX 

Starting trigger for the random number 
generator 

33333333333 

RNDOra 

Used to. adaptively alter (increase or 
decrease) the step size. for the random 
ray search 

.25 

INGRAN 

Number of sequential two sided unsuc- 
sessful random steps required before the 
step size will be reduced 

(AS = AS * RUDOWN) 

6 

lUPRAN 

Number of sequential two-sided success- 
ful random steps required before the step 
size will be increased 

(AS = (AS/RUD0i\71+ 2)/2.0) 

6 

RAN GEN 

Indicator for random number generator: 

>0, generates uniformly distributed random 
numbers 

.<0, generates normally distributed random , 
numbers . 

1.0 

FNORMAL 

Value by which the normally distributed random 
numbers will be divided in order to reduce 
the 3a value of the range over which random 
numbers will be generated. If FNORMVL=1.0, 
the 3a value will be approximately i tt 

5.0 
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NWEMONIC 

DESCRIPTION 

NOMINAL ' 
VALUES 

RBAR 

The R*s are the best historical estimate 
of the parameter perturbations required to 
produce a gain. The R’s are adaptively 
determined by the program and may be input 
for the purpose of ^^picking up** an optimi- 
zation run. 

100*0.0 

OLDRWT 

Weighting constant applied to the old R's 
when updating the R's following a gain 

5.0 

OLDGWT 

Weighting constant applied to d<()/ds when 
updating dcj)/ds following a gain. Where 
d(j)/ds represents the historical estimate of 
the gradient 

5.0 

RRSEXP 

Exponential constant applied to the ratio 
of current estimate of the gradient to the 
historical estimate of the gradient. The 
effect of this variable is to act as a 
dynamic weighting constant on the current 
random perturbation that produced a gain. 
The larger the gain the greater the effect 
on the R * s 

1.0 

RLIVT 

Weighting constant applied to the R*s in 
computing the current parameter pertur- 
bations . 

2.5 

I PARRS 

Directed random ray search parameter. 

0 

IPDRRS 

detail, and function values print indi- 

0 

IPFRRS 

cators 

0 


B.3.6.11 Arbitrary Ray Search (METHOP^ = 11) . The arbitrary ray 
search searches the ray passing tlirough two specified points in the 
multidimensional control space. For printout, it is suggested that 
the tabular summary feature of A.ESOP be used. 
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MNEMONIC 


DESCRIPTION 


NOMINAL 

VALUES 


RALOHI 

Defines which direction to search; i.e., 
TRUE, search "LO" "HI" 

FALSE, search "HI" ->■ "LO" 

.TRUE. 

NPTSRY 

Total number of evaluations to be used 
on ray search 

10 

RAYDIV 

Defines step size for ray search; i.e. 
control parameter step size = (XTENHIj^ - 
XTENLOi) /RAYDIV 

10.0,. 

XTENHIi 

Defines the control parameter values at 
the "HI"t end of the multidimensional ray 
to be searched 

100*0.0 

XTENLOj 

Defines the control parameter values at 
the "L0"t end of the multidimensional ray 
to be searched. 

100*0.0 


+ XTENHIj^ need not be greater than XTENLOi, and XTENLOj need not be 
less than XTENHIj 'in this search. These two arrays merely define 
the end points of a ray in the multidimensional control space. 


B.3.6.12 Jacobson Search (METHOPj = 12) . The Jacobson search algorithm 
is based upon homogeneous function rather than quadratic models. A 
consequence of this is that (N+2) step convergence is obtained for 
homogeneous functions of any degree and that no one -dimensional search 
is required. 


MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

JACMAX 

Number of Jacobson searches carried out 
when the Jacobson search is called. (Number 

50 


of Jacobson descents) 


ANGRAD 

If ANGRAD is .TRUE, the gradient vector . 
(3<J>/«ct.) is computed analytically by the 
system^model . If NAGRAD is .FALSE, the 
gradient vector is computed num.erically by 
AESOP using two-sided parameter pertur- 
bations 

.FALSE. 
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MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

PFACTR 

Controls the size of the parameter pertur- 
bations used for computing numerical 
gradient vector 

l.OE-7 


Aoi = PFACTR * - «lo^) 


IBX 

An array describing the parameter bounds 

=1, a^. is lower and upper bounded 

=2, is lower bounded only 

=3, oi is upper| bounded only 

-4, oi is unbounded -- 

50*1 

IPDJAC 

1 

Detail print indicator for Jacobson search 
=1, to generate detail print 

0 


B.3.6.13 Contour Mapping Search (NETHOP^ = 13) . The contour mapping 
option provides a two-dimensional graphic display of the performance 
function., The plane is selected by choosing any two parameters ai and 
aj . Plots are obtained on CALCOMP or Houston plotting devices. 


MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

IPDMAP 

Detail print indicator for mapper (Search 13) 

0 

NXMESH 

Number of x mesh points for map (maximum of 

0 


52 

0 

NYMESH 

Number of y mesh points for map (maximum of 
52 

0 

NZCUTS 

Number of contours for map 

b 

MAPALX 

Subscript of parameter for x mesh 

0 

MAPALY 

Subscript of parameter for y mesh 

0 

ZLMAX 

Maximum z level for contour map If ZLMAX= 
ZLMIN, then the program will use the max 
and min z mesh values 

0 

ZLMIN 

Minimum z level for contour map. If ZLMAX= 
ZLMIN, then the program will use the max 
and min z mesh values 

0 
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MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

LUMAP 

Logical unit to be used as the map plot 
output device. LUMA.P ^ Oj no plot will be 
produced; if LUMAP > 0 the program will 
read the following six cards for the map 
titles : 

0 


NCXT XTITLE > 

NCYT YTITLE ' ' ''' 

NCTl TITLEl y Format(3X,l2,lX,4A10) 

NCT2 TITLE2 j 

NCT3 TITLE 3 

NCT4 TITLE4 / 

/ /! 

f ^ " I . j 


Where NCXT, NCYT, NCTI,... etc. are the 
number of BCD characters in each title and 
XTITLE, YTITLE,... etc. are the arrays that 
will receive the BCD characters of the , 
titles. Each title may have a maximum of 
40 BCD characters. Figure 8 presents a 
contour map produced by AESOP and shows the 
placement of the titles on the plot. 

1 ' ' _ , ' 

; ■ \ 

ALTMAP. 

Logical variable used to indicate which 
bounds apply to the parameters'. If ALTM.AP 
is .TRUE, the bounds on the parameters . 

will be taken as XTENLCi- and XTENHIi. 
ALTMAP is .FALSE, the bounds on the 
parameters ui will be taken as ALPLOi and 
ALPHIi- 

.FALSE. 

XTENHIi 

Defines the upper bounds on the parameters 
when ALTMAP is .TRUE. 

100*0.0 

XTENLCi 

Defines the lower bounds on the parameters 
ai when ALTMAP is .TRUE. 

100*0.0 
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B.3.7 AESOP Print Control 


Program AESOP has a flexible print but capability. Varying levels 
of print out are available at the user's option. The print control for 
the individual searches was described in the previous section. This section 
will describe the general print controls that apply to all search algorithms. 

In all cases a value of one for the print control variable v;ill cause 
print out of the required data; a zero value will suppress the print out; 

B. 3.7.1 Input Data Print Out . The data input to AESOP via FORTRAN 
namelist lAESOP is out by a simulated namelist AESOUT. The output is 
organized by search algorithms and only printed for theuse'archeS. se.ljected 

by the user. • ■ . iroj iarrui- 


MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

IPNAML 

Controls print out of data input via FORTRAN 

1 


namelist lAESOP 


B.3.7. 2 

Print out of Function and Control Parameter Values 

During 

the Optimization Process. 


MNEMONIC 

DESCRIPTION ! 

NOMir^AL 



VALUES, 

IPINEV 

Will cause print out of the nominal values of 
the functions and control parameters for the 
initial report 

1 

IPDPEN 

Will cause print out of the function and control 
parameter values after each evaluation 

1 

IPGAIN 

Will cause print out of the function and control 
parameter values after each evaluation in which 
there was a gain in performance 

0 

IPSRCH 

Will cause print out of the function and control 
parameter values at the end of each search 

0 

IPACYC 

Will cause print out of the control parameter 
values at the end of each optimization cycle 

0 . 
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MNEMONIC 

. DESCRIPTION 

IPFCYC 

Will cause print out of the function values 
at the end of each optimization cycle 

IPFNEV 

Will cause print out of the function and 
control parameter values at the end of the 
optimization process (Final Report). 


NOi^JINAL 

VALUES 


3. 7':3o'i Tabular Sufamarvi Print ..Control . Tabular listing of the 
function and control pararaeter values can be produced to give the 
user a concise picture of the convergence history. 



DESCRIPTION 

NOMINAL 

VALUES 

LUCYC 

Tape or disk unit used to store the conver- , 
gence histories of the function and control 
parameter values for tabular listing at the 
end of each optimization cycle (cycle 
summary report). If LUCYC >< 0, no tabular 
listing at the end of a cycle will be 
produced. 

0 

LUFIN • 

Tape or disk \mit used to store the conver- 
gence histories of the function and control 
parameter values for tabular listing at the 
end of the optimization process (Final Report) 

0 

ITCYC 

1 

Controls content of cycle summary report: 

=0, the function and parameter values for 
all evaluations will appear in the 
cycle summary report ! 

=1, the function and parameter values 
will appear only for evaluations in 
which there was a gain in performance 
=2, the function and parameter values at 
the end of each search comprise the 
cycle summary report 

1 

ITFIN 

Controls concent of final summary report: 
=0, the function and parameter values for 
all evaluations will appear in the 
final summar)' report 

2 
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MNEMONIC 


DESCRIPTION 


NOMINAL 

VALUES 




=1, the function and parameter values will 
appear only for evaluations in which 
there was a gain in performance 
=2, the function and parameter values at 
the end of each search comprise the 
final summary report 

=3, the function and parameter values at 
the end of each optimization cycle' 

will comprise the final .summary report, 

• ' SJ. o 

; Y5iW3M 

LUCHDP 

Logical linit' (tape ,or disk) oh which the 
final summary report data is. to be refor- 
matted into arrays for plotting of conver- 
gence histories. If LUCHDP •< 0, the data 
is not saved 

0 


B3.8 Other AESOP Data 


This section describes the AESOP input data not defined in the 
previous sections. 


MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

IRLIST 

Controls the reading of the AESOP name- 
list data lAESOP. If IRLIST=0, no name- 
list data is input. This is used to 
suppress the input when all data is speci- 
fied internally by user's program 

1 

OPTIMA 

Not used by AESOP 

2 

DIALOG 

Logical variable used by AESOP. If 
.TRUE, the entire common for AESOP is to 
be saved and restored each evaluation 
by the main program 

.FALSE. 

EOF 

Logical variable. If .TRUE. AESOP will 
call EXIT and terminate the computer run. 

. FALSE . 
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MNEMONIC 

DESCRIPTION 

NOMINAL 

VALUES 

JJJ 

Evaluation counter 

-1 

MENTRY 

jTOCf^T 

Controls the dptimizatibh process. MENTRY is 

0 


set to nine (9) by AESOP when number of cycles 



exceeds TREPET or JJJ ^ MAXJJJ 
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PART II - PROGRAM STRUCTURE. TECHNIQUES, AND ANALYSIS 


5. COMBUSTOR PROGRAM STRUCTURE 


The general structure of the AUTOCOM program is designed around the 
physical and operational characteristics of the combustor, the weighting of 
these characteristics, and the rating of the system which is also based on 
weighting factor constants. Accordingly, the description for the following 
elements 

Combustor design inputs 
Specific combusfor’^ehafactefistics 
Average combustor characteristics 
Rating. equation 

is given in Section 5.1 through 5.4. 

Figure 5-1 depicts the task structure of the combustor program. As 
can be seen, the first set of designer inputs consisting of variable and 
fixed inputs describing the operational engine, is used to compute the 
specific combustor characteristics Fll through F82. Thereafter, the results 
of these 'Computations are weighted by a second set of designer inputs, 
neimely the weighting factor constants AFll through AF82, and BF72. These 
weighting factor constants reflect the designer's view regarding the im- 
portance of each individual specific combustor characteristic. Finally, 
the payoff criterion which is represented by the rating equation, combines 
the average combustor characteristics (FI through FVIII), by using designer 
inputs for the weighting factor constants AFI through AFVIII, BFI through 
BFVIII,and CFII. This class of weighting factor constants reflects the 
designer's view regarding the importance of each group of average combustor 
characteristics such as those for performance, stability, pressure drop, 
injector complexity, etc. The derivation of the weighting factor constants 
is. discussed in Section 5.5. 

Sections 5.1 through 5.4 present a discussion of the various elements of 
the program structure as shown in Figure 5-1. 


5.1 Design Inputs 

5.1.1 Combustor Design Variables (Input data for namelist COMDPV) . 
The independent design variables of the combustor are as follows with 
limits as indicated (for details see Section 2.3); 
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SPECIFIC COMBUSTOR 
CHARACTERISTICS • 


Fll THROUGH F82 

i5:.-:rBxO ? 


AVERAGE COMBUSTOR 
CHARACTERISTICS 


FI THROUGH FVIII. 


COMBUSTOR RATING 




OUTPUT 


Figure 5-1 

Task Structure of the Combustor 
Program 
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PROGRAM 

SYMBOL 

DESCRIPTION 

MAXIMUM 

VALUE 

MINIMUM 

VALUE 

BFLCIR 

Number of circumferential baffles 

Positive 

Integer 

0 

SBFLRAD 

Number of radial baffles 

.11 

0 

DC 

* Chamber, diameter, in. 

Positive 

Number 

Throat 

Diameter 

DF 

Diameter of fuel orifices, in. 

1 .T ' ■ » 

ff 

.01 

X 

o 

Diameter of oxidizer orifices, in 

i'i ^ 

tr 

.01. 


''"‘f tehg th^df ^com^ t ion chamber , . in. 

II 

.oi: 

LP 

^Length of fuel orifices, in. 

II* 

.01 

LOX 

Length of oxidizer orifices, in. 

II 

.01 

NP 

Number of fuel orifices 

Positive 

Integer 

1 

NOX - 

Number of. oxidizer orifices 

II 

1 . 

OF 

Propellant mixture ratio, /fuel 

BF71 


VOLF 

Volume of fuel manifold, in^ 

; 

Positive 
Number ‘ 

0 " 

VOLX 

Volume of oxidizer manifold, in*^ 

1 

. 0 


f 


5.1.2 ~"' Fixed Coiiibus tor ''Design Parameters (Input, data for namelist COMMIS). 
The fixed combustor; design requirements are items whidi are nqt changed by 

the designer' for a given engine. This includes all items that* are required 
in those equations that could be different for combustors with different 
propellants or requirements (for details see Section 2.4). 

5.1.3 Weighting Factor Constants (Input data for namelist COMWGT). 

The "rating" of a combustor design is based on using a weighted average of 
all the c}>aracteristics of a given combustor. The weighting factors are 
constants used in equations to obtain average characteristics or a rating. 

The constants are introduced to allow the, designer tp incorporate his views 
regarding the imporliance or validity of one technique for obtaining a given 
characteristic vs. another technique for obtaining the same characteristic. 

For example, if the designer believed that only one technique was valid 
for predicting the performance of a given design he would assign a unity 
value to the weighting factor constants for that specific characteristic 
and the constants for all the other performance characteristics would be 
zero. Similarly, the weighting factor constants in the equation to obtain 
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a single rating for a given combustor design are introduced to allow the 
designer to introduce the relative importance of different types of 
characteristics (for example, stability vs. performance). The weighting 
factor constants, therefore, gives the designer the same control and 
flexibility in the computer program as he has with current ”cut-and-try” 
design techniques. To establish a base point for the rating system a 
hypothetical ideal combustor is given a rating of zero. The hypothetical 
ideal combustor would have 100 percent of theoretical C* performance, 
infinite damping rate for all modes of instability, zero' chamber length, 
chamber diameter equal to the throat diameter, etc. The weighting factors 
are terms beginning with the letters A, B, or C followed by the symbol 
for the characteristic with which it is to be used. A list of the weighting 
factor constants is given in Section 2; 2, and their derivation for rating 
an engine is discussed in Section 5.5. 


5.2 Specific Combustor Characteristic's 

The specific combustor characteristics are calculated from well-defined 
equations or curves that have been published and accepted by the professional 
community. The specific combustor characteristics are denoted by a capital 
F to designate a characteristic, followed by a two digit number. The first 
digit corresponds to the Roman Numeral which designates the average engine 
characteristic in which the specific characteristic will be used. The 
second digit then designates which characteristic is being considered. The 
specific combustor characteristic is calculated from an equation containing 
combustor design variables and fixed combustor design requirements. The 
following specific characteristics are available in the program code: 

Fll is the percent mass vaporized of fuel 

F12 is the percent mass vaporized of oxidizer 

F13 is a C*-efficiency based on the mixing model 

F14 I Available for program extension. The former performance 
^ computations based on A. D. Little Correlations have been 

Fisj eliminated. They are, however, still included in Sections 
4 and 5.5 

F20 is the chugging decay rate based on the fuel system 

F21 is the chugging decay rate based on the oxidizer system 

F22 is the stability characteristic for pulsed operation based on a 
correlation 

F23 is the stability characteristic for npn-pulsed operation based on 
■ a correlation c 

1 

F24 is a stability decay rate characteristic for fuel 

F25 is a stability decay rate characteristic for oxidizer 
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F26 is the stability decay rate characteristic based on the sensitive 
time lag model for a longitudinal mode 

F27-is the stability rate characteristic based on the sensitive time 
lag model for transverse modes 

F28 is the stability decay rate characteristic based on the response 
function approach 

;F29 is the stability characteristic based on the non-linear stability 
analysis and includes various modes of oscillation • 

F31 is the fuel pressure drop characteristic • 

F32 is the oxidizer pressure drop characteristic 

F41 is the number of : fuel plus oxidizer holes characteristic 

> F42 is the .volume of the oxidizer dome characteristic 

F43 .is the volume of the fuel dome characteristic 

F44 is the. length of the oxidizer holes characteristic 

F45 is the, length of the fuel holes characteristic ^ ^ 

F46 is an injector type complexity characteristic 

F51 is the length to nozzel throat diameter ratio characteristic 

F61 is the chamber diameter to nozzel throat diameter ratio characteristic 

F71 is the mixture ratio characteristic 

F81 is the radial baffles characteristic 

F82 is the circumferential baffles characteristic 


The specific combustor characteristics Fll through F82 are computed in 
Subroutines SECFll through SECF82 respectively. These subroutines are 
discussed in detail in Section 7. 


5.3 Average Combustor Characteristics 

Average combustor characteristics (FI through FVIII) are used in the 
rating equation (Section 5.4) togive the optimum combustor design; they 
are computed in Subroutines ACFI through ACFVIII. Average characteristics 
are specified by a capital F to denote a characteristic, followed by a 
Roman Numeral to designate which characteristic. The average characteristic 
is a weighted average of several specific combustor characteristics calcu- 
lated by different techniques.. 
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5.3.1 Average Performance Characteristic. FI is the average performance 
characteristic based on the vaporization efficiency and the c -efficiency, and 
varies from 0 to 100 percent. The average characteristic’ equation is 

FI = AFll (100-Fll) 

+ AF12 (100-F12) 

+ AF13 (100-F13) 

+ AF14 (100-F14) 

+ AF15 (100-F15) (5-1) 

■ 5.3.2 Average Stability Characteristic. FII is ’• -tW •ivbr'a’geiSral^i r^tys’- ' i- 
characteristic based on an equivalent damping rate- and raries from' - 
(damps at an infinite rate with time) to + °° (grows at an infinite- rate 
with time) . The average stability characteristics equation is 

h=29 

e PII = ^5 e AF(n)*F(n) 

■ n=20 ’ : y (5-2) - 

or 

FII = ALOG[EXP(AF20*F20) 

+ EXP(AF21*F21) 

+ EXP(AF22*F22) 

+ EXP(AF23*F23) ... 

+ EXP(AF24*F24) 

+ EXP(AF25*F25) 

EXP(AF26*F26) 

+ EXP(AF27*F27) 

+ EXP(AF28*F28) 

EXP(AF29*F29)] (5-3) 


5.3.3 Average Pressure Drop Characteristic . Fill is the average 
pressure drop characteristic based on the pressure drop across the injector 
face and ranges from 0 to a positive number. The average stability 
characteristics equation is 


AF31*F31 + AF32*F32 

Fill = (5_4) 

AF31 + AF32 
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5.3.4 Average Injector Complexity Characteristic . FIV is the average 
injector complexity characteristic based on the number of injector elements, 
type of elements, injector dome volume and injector face thickness, and 
varies from 0 to a positive number. The average injector complexity character- 
istic equation is 


FIV = AF41*F41 + AF42*F42 

+ AF43*F43 + AF44*F44 

+ AF45*F45 + AF46*F46 (5-5) 


5.3.5 Average Chamber Length Characteristic . FV is the average 
chamber, dength, .characteristic based on the chamber length to nozzel throat 
diameter ratio, and varies from 0 to a positive number. The average chamber . 
length characteristics equation is = 

FV = AF51*(F51**AF52) (5-6) 


5.3.6 Average Chamber Diameter Characteristic. FVI is the average 
chamber diameter characteristic based on th¥ chamber diameter at the injector 
face to nozzle throat diameter ratio, and varies from 0 to a positive number. 
The average chamber diameter characteristics equation is 

FVI = AF61*[(F61-1)**AF62)] (5-7) 


5.3.7 Average Mixture Ratio Characteristic . FVI I is the average 
mixture ratio characteristic which varies from 6 to a positive number. The 
average mixture ratio characteristics equation is 

FVI I = AF71*[(BF71-F71)**AF72] (5-8) 

5.3.8 Average Baffle Characteristic. FVIII is the average baffle 
characteristic which varies from 0 to a positive number, and is 

FVIII = 2*L50(AF81*F81- +: AF82*F82) . * ..(5-9) 

where L50 = length of combustion chamber 
required to vaporize 50^ per 
cent of the propellant. 

If AF81 and AF82 = 0, the computation for FVIII (Subroutine ACFVIII) will 
be bypassed. 
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5.4 Rating Equation 


The specific equation to obtain a single numerical value, 0, for a 
combustor rating is given by 

<P = AFI(FI)®^^ 

+ AFlI*e(^PII-"CFn*FII) 

+ AFIII*FIIlC®Pi^^^5 
+ AFIV*FIV^®^^'^i 

+ AFV*FV*^®^^^) . . cc ifc Jisjasr 

+ AFVI*FVr^®^^^^ bnc-'o-r. ^ 7^3 '.tov 

+ AFVII*FVII^®^'^^^^ 

+ AFVIII*FVIII^®^^^^^^ C5-10) 


5.5 The Derivation o£ the Weighting Factor Constants for 

Rating an Engine 

The combustor rating is based on the payload penalty of an envisioned 
rocket using this engine. Payload penalties mean an increase in rating. 

The derivation of the weighting factor constants takes the "Engine Test Case" 
of Section 4 into consideration. However, this approach can also be applied 
correspondingly for other engines. Based on the 15000 lbs. thrust engine 
using H 2 /O 2 propellants, the following assumptions can be made: 

’ 1. Rating is in terms of payload. 

2. Cost of 1 lb. of payload is $4,000. (This is obtained by 
dividing total vehicle cost by number of pounds of payload.) 

3. Reducing cost of engine by $4,000 is equivalent to adding 1 lb. 
of payload to the rating. 

4. Adding 1 lb. of weight to engine is equivalent to taking a pound 
of payload out of the rating. 

5. Increasing the engine performance (specific impulse) by 1 second 
or 0.22 units of efficiency is equivalent to adding 33 lbs. of 
payload. (This is obtained by performing a mission analysis with 
various impulse values.) 

For the first term in the rating equation, based on the 15000 lbs. thrust 
engine: 

AFll is 0 because fuel is gas and already vaporized. 
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AF12 and AF13 are considered to be of equal importance (% LOX vaporized 
+ % mixed); therefore are assigned a value of 1. 

AF14 and AF15 are much less than AFll and AF12 beacuse the A. D. Little 
correlations are not considered to have a very high level of confidence 
in the prediction of performance. A 99% vaporized propellant can be 
considered equivalent to about a 1% value for F14 and F15. Therefore, 

AF14 and AF15 are given values of .01. 

BFI is unity because a trade study of payload with impulse shows a 
linear dependency between payload and impulse. 

AFI is 150 under the assumption of item 5 that for increasing performance, 
.22 units result in 33 pounds of payload. 

For the second term in the rating equation: 

F24 is not applicable to this engine as fuel is vaporized and, hence, 
Dykema’s vaporization model is not valid. Therefore, AF24 =0. 


All other AF20 to AF29 stability indices are considered of equal importance 
and significance; therefore, AF20, AF21, AF22, AF23, AF25, AF26, AF27, AF28 
and AF29 are assigned unity values. 

The AFII, BFII, and CFII weighting factor constants are obtained by the 
following three data points: 

FII = negative decay rate, then rating is 0. . 

FII = 0 (Neutral stability) is equivalent in importance to losing 

six seconds of impulse of 200 lbs. payload under the assumption 
of item 5, so rating is +200. 

FII = .10 (very unstable engine) is equivalent to 16 seconds of impulse 
or 527 lbs. payload, so rating is +527. 

For the third term in the rating equation: 

Analysis of turbopump has indicated that an increase in pressure or 
pressure drop of 200 psi will increase the weight of pump 8 lbs. and 
cost $8,000 for an equivalent rating of 10. An increase of pressure 
of 400 psi increases the weight 30 lbs. and cost $20,000 for an 
equivalent rating of 35.. Therefore, with Fill in units of psi: 


Fill = 0 

Rating = 0 

Fill = 200 

Rating = 10 

Fill = 400 

Rating = 35 
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From these conditions 

AFIII = .00069 
BFIII =1.8 

Pressure increase on fuel side results in a weight and cost penalty 
twice as large as LOX pressure due to low density of hydrogen. Hence, 

AF31 = 2.0 

AF32 = 1.0 

The fuel pump etc. for H 2 -O 2 system is approximated twice as big as 
the LOX pump, etc. Therefore, . it costs about; twice 'as'muchralsoT 

For the fourth term in the rating equation: 


For F41 it is estimated that increasing the number of holes from 
432 to 864 would increase the cost $12, -250 or an equivalent of 
3.06 lb. payload. Also it would increase the weight 1 lb. Therefore, 


AF41 = 3.06 +1 = .009 

432 

For F42 it is estimated that increasing Of the LOX volume by 
10 inches^ would increase the cost $1,000 and the weight by 1 lb. 

AF42 = .25 + 1 = . 125 

10 • 

For F43 apply the same results as for F42. 

For F44, a length increase of .4 inches is estimated to cost $1,000 
and increases the weight by 2 lbs. Therefore, 

AF44 = .25+ 2 =5.625 

.4 

For F45 a length increase of .06 inches is extimated to cost $200 
and increase the weight by .2 lbs. Therefore, 

AF45 = .05 + .2 =4.16 

■ .06 

For F46 the element type is not varied, so AF46 = 0. 
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'With respect to the rating equation it is reasoned that the rating goes 
as the square of FIV rather than the linear assumption used for F41,etc. in a 
region of doubling the size. Therefore BFIV = 2, 

Thus it follows, with FIV = 5^AF,^*F ^ ; 

^ (n) (n) 


FIV = 0 Rating = 0 


FIV = 8.625 Rating = X 

FIV = 17.25 Rating = X + 8.625 

For an engine with" FTVv= 17.25 by doubling engine size, rating increase 
is 8.625. Then 

X = 2.942 
AFIV = .0395 
BFIV = 2 

Since rating X = AFIV (8.625)^, it becomes 
X + 8.625 = AFIV(I7.25)^ 


For the fifth term in the rating equation: 

Under the assumption that FV is linear with F5I, AF5I and AF52 
(AF5I and AF52 = I) , all rating is taken up in AFV and BFV. 

To calculate FV coefficients the following data were used: 


F51 

Wt(lbs) . 

Cost 

Igp Loss 

Total 

0 

0 

0 

0 


1 

5 

10,000 

.1 . 

40.5 

2 

10 

■ 20,000 

4 

147 


Therefore BFV = 1.86 

AFV =40.5 


■« 
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For the sixth term in the rating equation: 

To calculate FVI coefficients, it is again assumed that AF61 and AF62 
are unity and all rating is taken up in AFVI and BFVI . The following 
data were used: 



Therefore, BFVI = 2 

AFVI =35 

For the seventh term in the rating equation: 

To calculate I from mixture ratio, the following data is used: 



*) According to the expression for FVII (Equation (5-8)), 
the magnitude of F71 as a variable can never reach the 
value of 5.06. 

Therefore BF71 = 5.06, then AF72 = 2 (to make a parabola). 


Arbitrarily set AF71 = 4, then 
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M.R. 

FVII 

- - - 

Rating 

5.06 

0 

0 

4.56 § 5.56 

1 

66 

4.06 § 5.06 

4 

330 


Then BFVII =1.16 

AFVII = 66 


Based on the preceding assumptions, the weighting factor constants 
used for the sample case of 'the 15000-lbs.- thrust engine in Section 4, 
are listed as follows: 


AFl = 150.0 
AFII = 200.0 
AFIII = .00069 
AFIV = .0395 
AFV = 40.5 


Constant of the performance characteristic in 
the rating equation. 

Constant of the stability characteristic in 
the rating equation. 

Constant of the pressure drop characteristic 
in the rating equation. 

Constant of the injector complexity characteristic 
in the rating equation. 

Constant of the chamber length characteristic 


AFVl =35.0 


Constant of the chamber diameter characteristic 


AFVII = 66.0 


Constant of the mixture ration characteristic in 
the rating equation. 


BFI =1.0 Exponential on the performance characteristic 

in the rating equation. 

BFI I = 0.0 Exponential on the stability characteristic in 

the rating equation. 


BFIII = 1.8 
BFIV = 2.0 
BFV = 1.86 


Exponential on the pressure drop characteristic 
in the rating equation. 

Exponential on the injector complexity charac- 
teristic in the rating equation. 

Exponential on the chamber length characteristic. 
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BFVI =2.0 
BFVII = 1.16 

CFII = .1 


Exponential on the chamber diameter characteristic. 

Exponential on the mixture ratio characteristic 
in the rating equation. 

Exponential on the stability characteristic in 
the rating equation. 


AFll = 0.0 
AF12 = 1.0 
AF13 =1.0 
AF14 = .01 
AF15 = .01 


Constant in the performance characteristic equation 
Constant in the performance characteristic equation 
Constant in the performance .'characteristic equation 
Constant in -the performance 'Characteristic^ equation 
Constant in the performance characteristic equation 


AF20 =1.0 
AF21 = 1.0 
AF22 = 1.0 
AF23 = 0.0 
AF24 = 0.0 
AF25 =1.0 
AF26 =1.0 
AF27 = 1.0 
AF28 =1.0 
AF29 =1.0 


Constant in the combustor stability characteristic 
equation. 

Constant in the combustor stability characteristic 
equation. 

Constant in the combustor stability characteristic 
equation. 

Constant in the combustor stability characteristic 
equation. . 

Constant in the combustor stability characteristic 
equation. 

Constant in the combustor . stability characteristic 
equation. 

Constant in the combustor stability characteristic 
equation. 

Constant in the combustor stability characteristic 
equation. 

Constant in the combustor stability characteristic 
equation. ' 

Constant in the combustor stability characteristic 
equation. 
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AF31 = 2.0 Constant in the fuel pressure drop characteristic 

equation. . . 

AF32 = 1.0 Constant in the oxidizer pressure drop characteristic 

r equation.' . 


AF41,.^„ . 009 . Constant, for the injector, orifice number. 

AF42;o?t j . 125 Constant., for the \Oxidizer dome volume. 

AF43r=? ai^ai'ioJ^.'.^Goni'tant-^for.- the fuel dome volume. 
AF44^o^:5V6^25^1^-- ’ 'Con’staht— for- the length of the oxidizer orifices. 

AF45'^-==^*4V16 Constant for the length of the fuel orifices. 

AF46 - 0.0 Constant for the injector type complexity. 

AF51 fe 1.0 Constant for the chamber length characteristic. 

AF52 5= 1.0 Exponential on the chamber length characteristic. 


AF61 = 1.0.. ; -Constant for the chamber diameter characteristic. 

AF62 =1.0 Exponential on the chamber diameter characteristic. 

Constant for the mixture ratio characteristic'. 

Constant for the ‘mixture ratio characteristic. 

Exponential on the mixture ratio characteristic 
equation. 

. i ' ■ 

Constant for the radial haffle characteristic. 

Constant for the circumferential baffle characteristic. 


AF71 = 4.0 
BF71 = 5.06 
AF72 =2.0 

AF81 = 0 
AF82 = 0 
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6. DESCRIPTION OF AUTOCOM PROGRAMMING TECHNIQUES 


Implementation of the application of multivariable search techniques 
to optimum combustor design follows three steps which form the basis for 
the computer program design: ' ' ' 

1. Combustor Point Design . A modular approach' permitting user . 
understanding and subsequent refinement, and rapid extension of the point 
design computer code, is employed in conjunction with algebraic FORTRAN IV 
coding throughout the program. 

2. Combustor Optimization . Optimization is performed using''the existing 
multivariable search program AESOP (Automated '^ngind'erihg -and j^ci.eritific' . '''' 
Optimization Program) ■, Appendix B*) This ■program*'can“rap’i'5iy 'lie'xdupIed\''=--'^'''“ 
to any point design program where the performance criterion i.s- a computable •. 
function of a finite number of design ’controlled parameters here ;a cqmp^li- . ^ 
cated weighted function of a set of specified combustor' char acfbristics^’ * *'■' 

3. Payoff Criterion . Combustor design performance comparisons are 
based on a single numerical value combining a set of average engine charac- 
teristics in a weighted fashion; in this case, the rating equation. The 
weighting factors are set to nominal values within the program. The designer 
is able to modify the weighting factors through the data input. 

To implement the three steps, the program is, in principle, synthesized 
by using the following building blocks: 

* Input/Output Control Program 

• ■ t j 

> Combustor Characteristics Computation Program (including average" 

■ combustor characteristics) 

* Combustor Rating . . , * . 

* Optimization Program AESOP 

* General Subprograms 

All progra.m coding utilizes the algebraic FORTRAN IV language. FORTRAN 
subset common to the following major digital computer systems is employed; 

IBM 7040/7094 

IBM 360 Series 

UNIVAC 1107/1108 

CDC 6400/6600 


*) Appendix B is included in Part I 
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6.1. Modular Programming 

Modular programming is. applied. Essentially, it consists of the use 
.of a y/ell ordered tree of subroutines which utilize a system of user- 
(combustor designer) oriented names to perform user-oriented functions. 
Subroutines communicate with each other through a single **C0MM0N” block 
containing, all data necessary for evaluation of the combustor rating once 
a set of combustor design variables have. been* specified. 


6.1.1 Use of COMMON . The combustor design program module utilizes a 
.^single named "COMMON ..which contains all combustor data, specific engine 
''characteristics , ; weighting.;cpnstants , average engine characteristics, fixed , 
design ; variabiesf, free f,design .variables , and combustor rating.. These 
quantiti escarp pres^e^ed an a. single. ,array;a^cess to . any item of information 
in the^^arra^ technique. A typical example illustra- 
ting ^this\pracedure^ .... . 


SUBROUTINE SECFll 
COMMON/COMBST7C®ATA(200O) 
EQUIVALENCE (LC, CMDATA(123)) 


RETURN 

END 

Here, the subroutine "SECFll" requires as input the combustor length, LC. 

This information is contained in the combustor design program array "CMDATA" 
at the (123) rd position. This method of communication between subroutines 
is rigidly followed. A directory subroutine contains a complete listing 
of all data in the common block. Two hundred vacant positions are maintained 
at the rear of the common block for subsequent program expansion. 


6.1.2. Combustor Program Outline . The combustor main program controls 
the combustor analysis and optimization procedure through a sequence of major 
subroutines designated as subprograms. Each of these subprograms performs a 
gross function such as input, output, optimization, specific engine charac- 
teristic evaluation, etc. A schematic diagram of the combustor main program 
illustrating major subprograms and their function is presented in Figure 6-1. 
Each subprogram defined in this figure consists of a number of subroutines. 


6 . 2 Input/Output Formats 

6.2.1 Input Format. An input procedure employing user-oriented 
mnemonics is“"applied. This approach avoids the tedious, rigidly formatted 
data input procedure employed by many existing FORTRAN IV and symbolic coded 
programs. The input procedure employed is the FORTRAN IV "NAMELIST" option. 
There are nine places for the input data setup of which two are optional. 

For more details see Section 2. 
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Figure 6-1. Schematic of Combustor 
Optimization Procedure 
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6.2,2 Output Formats . Several output data formats are supplied. 

1. Data Output. Relevant data output are printed at the beginning 
of a case to provide a record of assumed combustor data. 

2. Design Iteration Output. A running summary of major combustor 
characteristics is provided in the design optimization loop to 
aid the designer in interpretation of combustor design progress, 

3. Final Design Output. A detailed description of final combustor 
design variables and characteristics following completion of the 
design iteration loop is available. 

4. Optional Data Retrieval Output. At the designer *s option all 
fixed combustor design data including tabular properties are 
displayed. This gives the ^designer the ability to verify all data 
bui^lt into' the combustor design program and, in combination with 

„ . ther “input '.package, the ability to modify the fixed data in a par- 
ticular study, ’ ' . 

For more details about program outputs see Section 3.'; 


6.3 Combustor Characteristics Computation Program 

The specific combustor, characteristics subprogram (CMBSEC) will be 
described to illustrate the programming philosophy. This subprogram computes 
the specific combustor characteristics Fll to F82. Each of these specific 
characteristics is computed in a separate subroutine. The subroutines 
are named by a six- letter word commencing with the letters SEC (Specific 
Engine Characteristic) and terminating with three characters indicating 
which specific engine characteristic is computed! Thus, for example, sub- 
routine SECFll will compute the per cent mass- vaporized of fuel. The 
structure of subprogram CMBSEC is shown schematically in Figure 6-2,. ' It 
consists of a sequence of calls to working level subroutines each of which 
computes a particular specific engine characteristic. These working level 
subroutines vary in complexity from simple equations to complex tabular 
data computations. Each working level subroutine contains logic to bypass 
the particular specific engine calculation when the designer has indicated 
this requirement through the input weighting constants (see Figure 2-3) . 

Thus, for example, if the weighting constant AFI is zero, subroutines 
SECFll to SECF13 will return without computing Fll to F13, thus saving 
unnecessary computation. Again, if the weighting constant AF13 is zero, 
subroutine SECF13 will independently return without computation of F13 since 
the designer has indicated this particular specific engine characteristic is 
to be omitted from the rating equation (see Appendix A). 

A similar structure to that described for specific engine characteristic 
computation is employed for average engine characteristic computation. 
Subprogram ’^CMBAVC’^ thus calls eight subroutines, ACFI, ACFII, . . , ACFVIII, 
respectively, to compute the eight average engine characteristics. Again, 
computation by-passes on appropriate weighting constant values are provided 
(see Appendix A) . 
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Subprogram CMBSEC 


ENTER 



— Per cent mass vaporized of fuel 
calculated via NASA TR-67 


— Per cent mass vaporized' of oxidizer 
calculated via NASA TR't;67.^^ ' - - 


— C* efficiency determined by the 
mixing model of TN-2881 


— Available for extension 


— Available for extension 


Chugging decay rate based on the 
fuel system 


— Circumferential baffle 
characteristic 


RETURN 


Figure 6-2. Schematic of Specific Engine 

Characteristic Subprogram "CMBSEC" 
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6 . 4 Combustor Rating 


The combustor rating as the payoff criterion provides a single numerical 
measure of the combustor* s capability and is constructed on the basis of a 
weighted sum of the average engine characteristics. The weighting factors 
employed in computing the combustor rating are user-defined in the AUTOCOM 
code (see Section 5.5). These weighting factors are based on the impact 
of each average engine characteristic on vehicle payload capability; they 
define the payload penalty associated with each characteristic. The 
combustor design optimization process is based on minimization of the 
combustor rating and, hence, the payload penalty. The rating is clearly a 
function of the combustor design variables, and the weighting factors 
entering into both the rating equation and the average engine characteristics. 
In the- present general AUTOCOM coding, the values for the weighting factor 
constants are setot6'"**P- of *-*0**. The combustor rating is computed in 
subroutine CMBRAT using the rating equation given in Section 5,4. 


6.5 Combustor Optimization Progr^ 

6.5.1 The General Optimal Design Problem for AUTOCOM . The optimal 
combustor design problem ia essentially a large scale, non-linear multi- 
variable optimization problem. The independent variables are the geometric 
and physical parameters defining the combustor in detail- -the combustor 
design variables. The dependent variables are the specific combustor 
characteristics --the combustion efficiencies, statbility characteristics, 
etc. The specific combustor characteristics in turn define a higher order 
set of dependent variables, the average engine characteristics, which 
define a payoff function through the rating equation. 

In the general design problem. Reference 6-1, a set of design 
variables, a, are defined. A unique set of design characteristics, F, 
are obtained from the vector of design variables 

F = F(a) (6-1) 

In the design process one of these characteristics is selected for 
minimization or maximization. This is the payoff function 

ct) = cj)(a) (6-2) 

In some designs the payoff criteria to employ will not be self-evident to 
the designer. In this case he may seek to define value function, V, 
which involves some combination of the performance characteristics 

V = V(F) (6-3) 

The value function is then employed as the payoff function. In the optimal 
combustor design problem, V is specified by the rating equation as a 
weighted combination of average engine characteristics. An alternate 
approach to the use of the rating_jequation is that of seeking constrained 
extremal. Constraint functions, are selected from the average engine 
performance characteristics. These constraints can always be defined 
such that 


149 


f = >Ka) = 0 


(6-4) 


With this approach, for example, the designer seeks the optimum combustor 
design on the basis of one set of average engine characteristics while con- 
straining other average engine characteristics to acceptable values. 

Both methods, the value function approach and that of the constrained 
extremal, are included in the combustor design procedure. A general tech- 
nique for incorporating constraints into the optimization formulation is 


the well-known "penalty 

function" approach. 

An augmented payoff function. 

4>, is constructed by 

• M 



- i = (j> Will;. " 



i=l . 



Here, the W. are a set of positive constraint weighting factors.’ Provided 
the are sufficiently large in magnitude, minimization of Equation (6-5) 
corresponds to minimization of Equation (6-2) in the presence of the constraint; 
of Equation C6- 4) . In practice, the Wi are automatically determined within ‘ 
the optimization algorithms employed in adaptive fashion on the basis of 
constraint behavior. 

The designer may wish to impose inequality constraints on the design. 

For example, he may seek the best combustor design, based on a weighted com- 
bination of efficiency and stability characteristics with constraints on 
pressure drop, injector complexity, chamber length characteristics, chamber 
diameter characteristics, and propellant mixture ratio characteristics. In- 
equality constraints can be imposed on the formulation by a transformation . 
into equality constraint form. Equation (6-4) . Suppose an inequality is to 
be placed on the i^^ performance function; then define, a constraint, 
such that 

'I'j = Fi" ; F. > 0 

= 0; F. < 0 (6-6) 

Constraining ipj to zero is now equivalent to the constraint F^ < 0 . 

Frequently the designer imposes inequality constraints directly on the 
design parameters. (Thus, for example, the fuel orifice diameter must be 
greater than .010 inches). These limits are dictated by a priori knowledge 
of the combustor and its operating environment. Generally, then, the design 
parameters are subject to lower and upper limiting values, ^ and a", such 
that 

< a < (6-7) 

These constraints limit the region of feasible designs to a hyper -rectangle 
lying in the multidimensional design parameter space. Equations (6-1) through 
(6-7) define in symbolic fashion the optimal combustor design problem. 
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Conceptually, they define most industrial design problems; for, in practice, 
a designer must always seek to express his problem in terms of a finite number 
of parameters, 

6,5,2 Optimal Design Techniques , Methods for solution of non-linear 
multivariable optimization problems have received considerable attention 
during the sixties. In general, solutions are obtained by the iterative 
search procedures which are collectively becoming known as ’^optimal seeking 
methods." The increasing interest in these techniques stems both from their 
ready application through the digital computer and the ease with which their 
theoretical basis can be grasped. Generally, -the non-linear optimal seeking 
method has its basis in logic rather than the higher branches of analytic 
mathematics. In essence, the technique corresponds quite closely to the tradi- 
tional design cycle. Parameters are perturbed; the system is evaluated, 
and, on the basis of resulting performance characteristics, a new design is 
evolved. Figure 6-3 . presents a schematic diagram, of the optimal seeking 
approach. A nominal design variable vector aQ, is supplied tp the optimization 
algorithm. The optimizer, in turn, supplies the design parameter values to a 
digital model of the system being designed--in the present study a digital 
combustor point design evaluation module. This system functions, in "black 
box" fashion and returns the corresponding performance characteristics, F, 
to the optimizer. Based on inspection of these characteristics, which define 
payoff (value function) and constraint values, a new design, a, is supplied 
to the system, and the process repeats in iterative fashion until the optimal 
performance (f> = (j)* for the constraint levels ijj = D is attained. 

It can be seen that the optimization process is largely divorced from 
the system model. This fact permits ‘construction of generalized optimization 
programs which can readily be coupled to digital system models. These models 
may be expressly constructed with this object in mind, or, equally, they may 
be existing digital system models constructed for conventional designer control 
and perturbation. An example of this type of generalized optimization is 
program AESOP, (Automated Engineering and Scientific Optimization Program), 
References 4-1 and 4-2 recently constructed under contract to the National 
Aeronautics and Space Administration’s Office of Advanced Research and Tech- 
nology. A detailed discussion of Program AESOP can be found in Appendix B, 

Part I . 

6.6 General Subprograms 

Evaluation of specific combustor characteristics requires the’ use of a 
variety of standard numerical subprograms for tabular data storage, table 
look-up, non-linear complex equation solution, etc. \The Aerophysics Research 
Corporation program library, which consists of approximately 500,000 cards 
containing all basic numerical subprograms required in construction of the 
combustor design procedure, has been employed for these basic operations . 

The four major subprograms whose positions in the master chart. Appendix A, 

Part I, are indicated and are listed in the following: 


151 




Figure 6-3. Optimizer Schematic 
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Subprogram Name 

Description 

CMBINI 

Initializes all fixed data in 
the combustor data base 

CMPRED 

Calculates and/or stores all 
fixed combustor design character- 
istics . 

FITCRV 

Fits curves for specific impulse, 
combustion temperature , specific 
heat ratio and molecular weight 
of the combustion gas as a function 
of the propellant mixture ratio 
and stores this information in the 
data base. 

POSTAL 

Computes minor design variables 
dependent functions . 


6.6.1 Subprogram CMBINI . This subprogram initializes all fixed data 
in the combustor data base. Several kinds of data are stored in subprogram 
CMBINI: . 

1. Constants such as: 

-Phase angle due to delay in burning in Subroutine SECF28, 

F28TBS =1.9 

-Acceleration of gravity^ in/sec^ , GREF = 32.174*12 

-Fuel injection distribution eccentricity in Subroutine 
SECF23, IDEF = 1.0 . 

-Log of fuel injection distribution in Subroutine SECF23, 

LIDF =0.0 

-Constant PI = 3.14159 

-Universal gas constant, in-lb/mole^R, RGUNIV = 1544*12 

“Vapor pressure liquid temperature coefficient in Subroutine ♦ 

SECF28, TEMCOE =8.0 

-Turbulence intensity in Subroutine SECF13, %, TURBIN =3.0 

-Specific heat ratio for hydrogen in Subroutine SECF28, 

SPHTH2 =1.2 

2, Physical constants for heptane as required for the scaling relation 
presented in Equation (7-2). The following table depicts the heptane 
values stored: 
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Program 

Symbol 

Value 

Description 

MUC 

1.74*10-5 

Viscosity of heptane, Ib/in.-sec 

RHOC 

2.47*10-2 

Density of heptane, Ib/in^ 

SIGLC 

1*10-^ 

Surface tension of liquid heptane, 
Ib/in. 

TLOC 

528.6 

Temperature of heptane, or 


3. The propellant interaction index, n, for non-hypergolic propellants 
with coaxial injection in Subroutine SECF26 

NTAUNl = 0.45 

and the propellant interaction index, n, for hypergolic propellants 
in Subroutine SECF26 ' 

NTAUN3 = 0.70 

4. Modes of oscillations in Subroutine SECF27 if no baffles are present; 


SNUH(l) =1.8412 
SNUH(2) = 3.0543 
SNUH(3) = 3.8317 
SNUH(4) = 7.0156 
SNUH(5) = 5.3313 
SNUH(6) = 8.5263 
SNUH(7) = 6.7060 

5, Argument of Bessel function in Subroutine SECF28 if no baffles are 
present: 


BESARG(l) = 1.8412 
BESARG(2) = 3.0543 
BESARG(3) =4.2012 
BESARG(4) = 3.8317 
BESARG(5) = 0.0 (longitudinal) 


6.6.2 Subprogram CMPRED . This subprogram calculates and/or stores 
fixed combustor design characteristics such as: 

1. Constant for number distribution in Subroutines SECFll and SECF12, 
ADIS =1.0 

2. Combustion chamber throat area: ATHRTI = (PI*DTHRTI**2)/4.0 
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3, A matrix of the relationship for the injector element type, ETYP; 
the propellant type indicator, PTYP; the propellant combination 
- parameters F02 and F03; and the principal element type PEI and PE2 
I required for Subroutines SECF22 and SECF23 

^4. The reduced initial propellant temperature for fuel, TLORF,and for 
oxidizer, TLORX, as used in Subroutines SECFll and SECF12: 

TLORF = TLOF/TCRF 


and 

TLORX = TLOX/TCRX 

(for definition of the program symbols see listing of 'Section 2.4) 


6.6.3 Subprogram FITCRV . Subprogram FITCRV is used for optimization runs 
to determine the. .four parameters: 

Isp = Specific impulse 
k = Specific heat ratio 
M = Molecular weight of combustion gas 
T^ = Temperature in combustion chamber 

as a function of the propellant mixture ratio 0/F. The procedure for storing 
new propellant data is described in Section 2.7.1. 

Subprogram FITCRV is a simple and rapid degital computer code used to 
find the polynomial Pj<^(x) of degree ^ k which best approximates, in the least 

squares sense, a weighted set of data points, x^ , y^ , with i = 1, 2, N* 

using orthogonal polynomials. Moreover, this subprogram calculates the 
coefficients of the fitted polynomial and the fitted values for a given set 
of arguments. . ^ 

For the propellant combination H 2 /O 2 as an example, the four propellant 
polynomial functions, k, M, and T^^, are nominally established within the 

program as the following^functions of the oxidizer/fuel ratio, 0/F: 

= 235.436554 + 121 . 376531*0F 

bp 

- 24.6650523*(OF)2 

+ 2. 09800922* (OF) ^ 

- ,0733124871*(OF)'^ 
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k = 1.46481804 - . 106859541*OF 
+ .0213594157*(OF)^ 

- .0021301452*(OF)^ 

+ (83.0205334*10-6) *C0F)4 

Tj = 3022.96538 - 554. 118155*OF 
+ 478. 420285* (OF) 2 

- 71. 4672888* (OF) ^ 

+ 3.19696006*(OF)4 

M = 1.39758043 + 3 .03152880*OF 

- .530108602*(OF)2 
+ .107300022* (OF) ^ 

- .00750532403*(OF)"^ 


When tabular data is entered, the subprogram will fit the polynomial curves 
through the tabulated points by a least squares procedure. Maximum poly- 
nomial order which can be specified is ten. In the case of H^/O^ propellant 
properties, this could be adequately represented by fourth order polynomials 
The ^goodness of fit" is illustrated by Figure 6-4. Here, the input curves 
are compared to the polynomial approximation at two sets of points. The 
circled points represent input tabular points; the crossed points represent 
evaluations midway between the input points. 

A broader treatment of the computer code for curve-fitting can be found 
in Reference 6-2. 


6.6.4 Subprogram POSTAL . Subprogram .POSTAL computes, characteristics 
and minor design variable dependent functions which are inside the optimiza- 
tion loop. A detailed discussion of the computations performed in Subprogram 
POSTAL is presented in Section 7.9, "Minor Design Variable Dependent Functions." 
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SPECIFIC IMPULSE. loD COMBUSTION TEMPERATURE,!, , (‘’f) 



PROPELLANT MIXTURE RATIO, o/F 


Figure 6-4. Polynomial Propellant 
Properties 




7. OUTLINE OF THE. AUTOCOM PROGRAM ANALYSIS 


The AUTOCOM program automatically determines the combustor characteristics 
given the chamber design variables and fixed input values. The analysis 
considers performance, stability, combustion chamber geometry , injector 
complexity and baffle characteristics. In an optional mode of operation, 
the program possesses the ability to automatically perturb the design para- 
meters defining the engine characteristics (optimization) . Performance and 
stability analysis modules (subroutines) as well as several other combustor 
characteristics modules (subroutines) available within the program are listed 
below and described in the following: 

. * Performance characteristic modules 

* Stability characteristic modules 

* Pressure drop characteristic modules 

* Combustor complexity characteristic modules 

* Combustor chamber length characteristic modules 

* Combustor chamber diameter characteristic modules 
Propellant mixture ratio characteristic modules 

* Baffle characteristic modules 


7.1 Performance Characteristic Modules (Subroutines) 

There are three subroutines, SECFll through SECF13, which compute perfor- 
mance related characteristics. In addition the two modules , SECF14 and SECF15, 
are available for program extension. * ! ~”- 


7.1.1 Subroutines SECFll and SECF12 . Subroutine SECFll computes the 
per cent mas^ vaporized of fuel, Fll, and Subroutine SECF12 computes the per 
cent mass vaporized of oxidizer, F12. If the propellant is in a gas phase 
when injected, the mass vaporized is assigned a value of 100 per cent (executed 
by input of ISTATF = 1, or ISTATX =1 in Input Namelist COMMIS) . Both 
subroutines are based on the same method using the same propellant vaporiza- 
tion subroutine VAPRIZ (Figure 7-4) and the polynomial curve fit routine 
YPOLYV. 

The values for Fll and F12 are determined from the technique described 
in Reference 7-1. The dxop size is first determined for heptane using 
Figure 7-1 (Figure 32 of Reference 7-1) . Options for parallel jets, 
impinging jets, and triplet jets are provided. The data of Figure 7-1 
indicates the drop size of heptane; for the three types can be determined 
from an equation involving the jet diameter, dj , of the form 
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Mass-median drop radius, 
i nches 


:0 



) - • ■ ' ■ . . 

Figure 7-1. Drop Sizes Determined from Experimental 
Engine Performance 
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(7-1) 


“ Adj® 

which generates a straight line on log-log scales. The constants A and B 
for the three types of jets shown are given in the table below. 



A 

B 

Parallel jets 

.0312 

.6194 

Impinging jets 

.046 

.9935 

Triplet jets 

.00575 

.6217 


Once the drop size for heptane is determined, the drop size for any other 
propellant, fuel or oxidizer, can be determined from the scaling relation ^ 


where 



(7-2) 


p = liquid density of the propellant at injector 
exit conditions, Ibm/in^ 


0£ = surface tension of the liquid propellant, Ibf/in. 

y = viscosity of the liquid propellant, Ibf/in.-sec 

The non-subscripted values represent heptane values while the subscript x 
denotes the input propellant values. From the drop size determined in this 
way and adding other propellant and chamber characteristics i’^a generalized' 
length can be calculated as shown below (Equation (103) of Reference 7-1) . 


gen 


where 


.44 


300 


. 66 ' 


Cl-Til,o,R) 


.4/ ""m 


.003 


1.45/ V ..75 /X N.8/M.at.35 

^ I-) I- 

1200 / 


140/ I 100, 


= chamber length , inches 


(7-3) 


A = chamber contraction ratio 
P- = chamber pressure ; Ib/in^ 


Vg = injection velocity^ inches/sec 
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^j^,o,R ~ initial propellant temperature ratio 

r^ = propellant drop size , inches 

X = heat of vaporization of liquid propellant , 
BTU/lbm 

M = molecular weight of liquid propellant, 
Ibm/lb.Mole 


a = constant for number distribution 


The generalized length approach was selected over the effective length 
approach as discussed in Reference 7-1 because of the applicability of the 
former approach to all propellants. With the generalized length calculation, 
the per cent mass, vaporized may be determined from Figure 7-2 (Figure 28f, 
Reference 7-1). This curve was fitted with a sixth order log-linear polynomial 
of the form 


PMV = Aq + log (£g^^) + A2[logC£g^^)]2 + ... + A^ [log(£g 3 „) 

where 


- 25.67 


-.394 

19.47 


-.007 

- 5.013 

^6 

.012 

= .381 




• ■ • . j ■ , ' . • 1 

The correlated results of this curve fit are also shown in Figure 7-2. The 
symbols denote values from the polynomial approximation. 


Figure 7-3 depicts the flow chart for controlling subroutines SECFll 
and SECF12. The computation- of subroutine VAPRIZ is presented in a more 
detailed flow chart. Figure 7-4. 



161 



162 


Figure 7-2. Correlated Results for Mass Vaporized is a Function 

of the Generalized Length (All Propellants, Operating 
Conditions, and Design Conditions, with a Geometric 
Standard Deviation of 2.3) 


YPOLYV. 


VAPRIZ 

POLYNOMIAL CURVE 
FIT SUBROUTINE 


EXECUTES COMPUTATION 
FOR PER CENT MASS 
VAPORIZED ^ 

(DETAILS IN FIG. 7-4) 



COMBUSTOR DESIGN 
LOOP 


SECFll 


COMPUTES PER CENT 
MASS VAPORIZED OF 
FUEL 


I 


SECF12 


COMPUTES PER CENT 
MASS VAPORIZED OF 
OXIDIZER 



CONTINUE COMBUSTOR 
DESIGN LOOP 


Figure 7-3. Flow Chart, for Subroutines SECFll and SECF12 
for Computing Per Cent Mass Vaporized of Fuel 
(Fll) and of Oxidizer (F12) . 
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ENTER 



RETURN 


Figure 7-4. Detailed Flow Chart for Subroutine VAPR I Z : 

(for Definition of Symbols see Equation (7-3)) 
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7.1.2 Subroutine SECF15 . Subroutine SECF13 computes the C*-efficiency in 
in per cent, F13, based on the mixing model as presented in Reference 7-2. 
Figure 7-5 (Figure 4 of Reference 7-2) depicts the relationship of the C*- 
efficiency vs. the mixing parameter, a . If a propellant combination is not 
in Figure 7-5, the curve for JP- 4 /O 2 can be taken. It is 

a = T*L^/S (7-4) 

where 


T = turbulence intensity 

= combustion chamber length 

S = average spacing between elements 

C* is found by a polynomial curve fit procedure YPOLYV so that F13 = 
YPOLYV ( ALPHA, PCCSTR, NPCSTR) . The polynomial coefficients PCCSTR and 
NPCSTR are stored in subroutine CMBINI . 

The flow chart for the controlling subroutine SECF13 is presented in 
Figure 7-6. 


7.2 Stability Characteristic Modules (Subroutines) 


Combustor stability characteristics are evaluated on the basis of the 
methods given in References 7-3 through 7-13. Several of these methods 
require the solution of * characteristic equations in the complex plane 
(Subroutines SECF20, 21, 26, 27, and 28). In certain cases, the neutral 
stability boundaries can be found in closed form. In the AUTOCOM method the 
actual decay rate rather than the neutral stability boundary is required; 
for the method must recognize how stable or unstable an engine is in order to 
define an improved design. 

A general technique for solving complex characteristic equations is 
embedded within the AUTOCOM code. This method is again based on multi- 
variable search, -this time in an inner optimization loop within the system 
evaluation module. 

Let a particular complex characteristic equation be represented by 

F. (?) = 0 

where ? is position in the complex plane, and let 

z = (aj, 

Then the roots of the characteristic equation can readily be found by solving 
the nonlinear optimization problem 
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C* efficiency, percent 


\- 

IV 



(b) Fluorine with various fuels. (c) Oxygen and fluorine with hydrogen. 



Mixing parameter, a 


Id) Oxygen and fluorine with ammonia. (el Hydrazine with various oxidants. 


Figure 7-5. Characteristic Exhaust Velocity Efficiency as a Function 
of Mixing Parameter at Oxidant-Fuel Weight Ratio for 
Maximum Theoretical Characteristic Exhaust Velocity 
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COMBUSTOR DESIGN 
LOOP 


YPOLYV 

POLYNOMIAL CURVE 
FIT SUBROUTINE 




COMPUTES THE C - 
EFFICIENCY AS DEFINED 
BY NASA TN-2881 


kA 


CONTINUE COMBUSTOR 
DESIGN LOOP 


Figure 7-6. Flow Chart for Subroutine SECF13 for 
Computing C*-Efficiency 
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$* = Miji[|F^(z) |] 

a . jj , ' 

"ii* 

Where multiple roots exist, each root obtained may be swept out of the > 
characteristic equation by multiplication by (z - Zj.qq^) > in this manner, all 
acceptable roots within a region of the con 5 >lex plane (the decay rate/ 
frequency- plane) may be found. 

The above method is particularly appropriate in a combustor design 
optimization in which a sequence of perturbed engine designs are evolved. 

For in every design evaluation following the nominal, the approximate loca- 
tion of stability roots are known from the previous engine design. This : , 
results in greatly improved convergence rates for location of perturbed 
stability roots on successive designs. 

Typical contours of the function |f.(^)| for the stability analysis 
are presented for illustrative purposes^in Figure 7-7. The effect of 
boundary selection in dicated. Thus, with the original search boundaries, 
shown dashed, no root is located; however, by increasing the search boun- 
daries the desired root is found. In evaluation of a nominal design, some 
manipulation of search boundaries is often required (for contour mapping - :r 
see Appendix B, Part I). 



Figure 7-7. Typical Roots from Stability Transverse Mode, 
Non-Coaxial Injection, Oxidizer 
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7.2.1 Subroutines SECF20 and SECF21 . Subroutine SECF20 computes the 
chugging decay rate, F20, based on the fuel feed system; and subroutine SECF21 
computes the chugging decay rate, F21, based on the oxidizer feed system. 

Both, F20 and F21, are determined from L. Crocco and S. I. Cheng’s theory 
as given in Refenence 7-3. 

For pump-fed systems. Equation (2.04.03) of Reference 7-3 is used in the 

form 


t- (1 + Es + JEs^)[l + s - n + ne"’’^] + PEse""^^ = F(s) = 0 (7-5) 

I 

This equation is coded in subroutine C20403. 

For pressure-fed systems. Equation (2.05.02) of Reference 7-3 is used in 
the form 

[1+Js + JEys^ + J^Ey(l-y)s (1+s-n + ne”^^) + Pe”^ (1+JEys^) = F(s) = 0 

(7-6) 

This equation is coded in subroutine C20502. 

The definition of the constants in Equations (7-5) and (7-6) is 


2APpx 


E = 

m 0 

g 

[ elasticity parameter ] 

(7-7) 

P 

■ 2AP 

[ pressure drop parameter ] 

(7-8) 

Z m 



j = 

[ inertia parameter ] 

(7-9) 


2APA0 

& 

The optimization program AESOP is used in an inner loop as a root finder 
by minimization of the function |F(s)| in Equation (7-5) or Equation (7-6). 
The Laplace transformation variable "s" has been solved for in terms of a 
real and an imaginary value. These are the elements of the complex root, 
defining decay rate and frequency of oscillation. 

The root searching procedure begins with activating subroutine STARTR 
which computes a neutral stability point for solving Equation (2.04.06) or 
Equation (2.05.04) of Reference 7-3 based on the coded-in guessed roots 
presented in Table 2-1. This is executed in subroutines C20406 and C20504, 
respectively. The neutral stability point thus found is used as the starting 
point in a search for the roots of Equation (7-5) or Equation (7-6). 

When a root is found, it is swept out by a root sweeping routine, SWEEPR, 
and another root solution is attempted. This process is continued for a pre- 
determined number of roots. The resulting complex roots are examined, and 
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the root with the largest positive real part and a non-negative imaginary 
part is selected. The characteristic function is the real part of the root 
devided by the gas residence time, 0 . Thus, 

s 


and 


F20 = Re(s) 


fuel 



(7-10) 


F21 = Re(s) /0„ (7-11) 

ox g 

where the imaginary part of the solution- represents the angular frequency. 

The symbols used in Equation (7-5) through Equation (7-11) are; 

, n = pressure index of interaction between combustion 
... pro.cesses and -oscillations in the combustion 

chamber 

' . ; i * • 

= combustion chamber pressure 

- . AP = injector pressure drop 

m = steady state mass flow rate 

^ = length of feed line (the orifice length of the 
injector for fuel or oxidizer is replacing 
"Jl" in this computer model) 

A = cross-sectional area of feed line (the total 
cross-sectional area for fuel or oxidizer holes 
■ is replacing. "A,', ' in this computer model) 

T = t/0 = sensitive time lag 

& 

^ = time lag : , 


For the case of supercritical or gaseous state 
of the fuel, the definition for t has been 
modified to the expression: 


, t 
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.G^ 


"^fuel 

and 


^SO(ox) 

Vs V. 

inj (fuel) 

Vo(fuel) 

7s V. , 

' 3 inj (ox) 


(7-l'2) 


(7-13) 


0g = gas residence time 
Cl) = angular frequency 


= length of combustion chamber required to 
vaporize 50 per cent of the propellant 

Vfnj = propellant injection velocity 

y = position of the equivalent concentrated line 
’ capacitance as fraction of feed line length. 
In the present program coding, y is 1,0 

X = compressibility factor of propellants 

p = propellant density 

s = Laplace transformation variable 


Note: if both, fuel and oxidizer, are in a gaseous state, 

the expressions for and x (Equations (7-12) 

and (7-13), respectively) are n8^ longer valid. 


The designer's choice for a pump-fed or pressure-fed system can be 
executed by input of the cards, listed below, in Data Input Namelist 
COMMIS (see Section 2,4): 


V ' 

' PUMP -FED 

PRESSURE -FED 

Fuel 

lEQF = 1, 

lEQF = 2 , 

Oxidizer 

lEQX = 1, 

lEQX = 2 , 


Figure 7-8 shows the flow chart for the controlling subroutines SECF20 
and SECF21. Figure 7-9 depicts the logic of subroutine STARTR, whereas 
Figure 7-10 shows a detailed flow chart for function LFIFTY. 
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Figure 7-8. Flow Chart for Subroutines 
. SECF20 and SECF21 for 
Computing Chugging Decay- 
Rates, F20(Fuel), and 
F21 (Oxidizer) 
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RETURN 


gure 7-9. Flow Chart for Subroutine STARTR to Comoute 
Neutral Stability. Starting Point for Solving 
Equations (2.04.03) and (2.05.02) in 
Subroutines SECF20 and SECF21 
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Figure 7-10. Detailed Flow Chart for Function LFIFTY which 
Computes the Length of the Combustion Chamber 
to Vaporize 50 per cent of . Fue.l or Oxidizer. 
(For Definition of Symbols see Section 7.1.1) 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF20 AND SECF21 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF20 and SECF21 
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7,2,2 Subroutines SECF22 and SECF23, Subroutine SECF22 computes the 
stability characteristics for pulsed operation based on the A, D, Little 
correlation (Table II of Reference 7-4) , F22, and subroutine SECF23 
determines the stability characteristics for non -pulsed operations, F23, 
also based on the A, D, Little correlation (Table I of Reference 7-4) , 

The correlation equations are shown in Tables 7-1 and 7-2, Here, SPIO 
indicates the occurrence of an induced, undamped, high or intermediate 
frequency oscillation, and SPIA indicates the occurrence of a spontaneous, 
high and intermediate frequency oscillation. 

The stability characteristics are obtained from: 


F22 = 200(SP10 - 

.5) 

(7-14) 

F23 = '200(SP1A - 

■.5) 

(7-15) 


SP 10.= [ - 5.55684 - 0.16675(LD - 1.45926) - 0.2657:(LDj - 8.20992) 

♦ 0.3456 (TPVM - 4.40439) ♦ 0.00009267(V^ - 1065.97697) 

- 13.678Sl(log D , - 1.69583) ♦ 0.15923(F02 - 0.12461) 

01 

- 74.95672(LR - 0.08022) - 0.36162(BF - 0.11682) 

+ 12.5715KLR - 0.08022) (PEI - 0.14330) 

♦ 0.27624(BF - 0.11682) (PE2 - 0.42991) 

- 0.003242(TPVM - 4.40439)^ > 0.42891(log D ^ - 1.69583)2 

of 

- 0.044068(LD - 1.45926) (TPVM - 4.40439) 

+ 0.0003075(V - i065. 97697) (log D . - 1.69583) 

o of 

- 3.0964KLR - 0.08022) (LDf - 8.20992) 

♦ 0.0003577(LR - 0.08022)(Vq - 1065.97697) 

- 163.91664(LR - 0.08022) (log - 1.69583) 

- 0.12184(PE1 - 0. 14330) (LD^ - 8.20992) 

. 0.15581(PE 2 - 0.42991) (LO - 1.45926) 

- 0.00002194(PE2 - 0.42991) (V^ - 1065.97697) 

+ 0.03835(PE2 - 0.42991) (TPVM - 4.40439) 

- 0.40008(PE2 - 0,42991) (log - 1.69583) ] 


TABLE 7-1. A. D. Little Correlation Equation to Determine 
the Stability Characteristics for Pulsed 
Operation 
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SPIA ~ [ 0.495S3 - 0.0014639(LR - 0.25065}(P^^ - 329.95025) 

♦ 0.024960(L - 16.24385) - 0.45096(iog R - 0.33429} 

I c 

♦ 7.15205(BF - 0.36833)^(log ID^ + 0.01452) 

+ 0.11250 (LR - 0.25068) QID - 2.13785) 

- 0.05177(F05 - 0.28507)(IDE^ - 1.92209) 

- 0.27789CF03 - 0.28507) 

- 0.0006173(IDE^ - 1.92209) (EPA - 3.76607) 

♦ 1.43317CLR - 0.25068)(MPE - 1.93154) 

♦ 2.09309CLR - 0.25068)’- 0.30534(BF - 0.36833) 

- 0.80153CPEL - 0.16561) + 0,001229(1^ - 16.24385)^ 

♦ 0.002199(MP0 - 1.93154)^ - 0.26849 (log R^ - 0.33429)^ 

4 0.00035321(P^. - 329.95025) (LD - 2.13785) 

- 0.b8700S(Lj - 16.24385) (log-R^ - 0.33429) 

- 0.17222(Lj - 16.24385)(log IDf + 0.01452) . • 

- 0.000131(Lj - 16.24385) (NIPE - 1.95154) 

- 0.81565(log IDj + 0.01452) (IDE^ - 1.92209) 

- 0.030990(LD- 2.13785)(MPE - 1.93154) ’ 

- 0.033444(IDE^ - 1,92209) (.*-IPE - 1.93154) . ' 

+ 0.25448(PE1 - 0.o6561)(LD - 2.13785) 

- 0.90924 (PEI - 0.16561) (MPE -“1.93154) 

0.17428(PE1 - 0.16561)(IDE^ - 1.92209) 

4 0.0957(PE1 -0.16561) (log R^ - 0.33429) 

- 0,000346S(F03 - 0.28507) (P^. - 329.95023) 

“ + 0.03245CF03 - 0.28507) (L^ - 16.24385) 

+ 0.0009S624(BF - 0.36833) (? . - 329.95023) *: 

Cl 

- 1.20213(LR - 0.25068)'(log R^. - 0.33429)) 




TABLE 7-2. A. D. Little Correlation Equation to Determine 
the Stability Characteristics for Non-Pulsed 
Operation 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF22 AND SECF23 
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7.2.3 Subroutines SECF24 and SECF25 . Subroutine S.ECF24 computes the 
decay rate for fuel, and Subroutine SECF25 computes the decay rate for 
oxidizer based on the high-frequency stability analysis of 0. W. Dykema 
(Reference 7-5] . He arrives at a dimensionless correlating parameter, N , 
called the stability number, which essentially represents the dimensionliss 
ratio of a characteristic molecular diffusion time to a characteristic 
acoustic time. The stability number is : ; 


where 


f P 

Ng=A[nd3][ —] (7-16) ■ 

w 

i 

A’ = stability constant J 

d = injector orifice diameter (fuel or oxidizer), in. 


f = frequency that gives maximum value for per cent 
combustion gain, cps 

n = number of injector orifices (fuel or oxidizer) 

= combustion chamber pressure at injector end, psia 
w = propellant flow rate (fuel or oxidizer) , Ibm/sec 


Because the constant A is very difficult to evaluate for, an operating thrust 
chamber, it must be determined by correlation of instability data. An 
average value of 6.0*10"^ has been assumed for this study. Thus, the equation 
for can be rewritten 


f 

N = 6.0*10"^*p(n*d^) (— : (7-17) 

^ w 


where pis the specific gravity of the propellant. 

The first bracket in the above equation includes the only terms which 
define the injector geometry. The second term includes the operating conditions 
of chamber mode frequencies, chamber pressure, and propellant flow rate. 

I- 

Frequency modes which have been considered are: 


* Frequencies of the 

* Frequencies of the 

* Frequencies of the 

* Frequencies of the 

* Frequencies of the 


1st, 2nd, and 3rd longitudinal mode 
1st, 2nd, and 3rd tangential mode 
1st and 2nd radial mode 

1st longitudinal with the 1st tangential mode 
1st longitudinal with the 2nd tangential mode 
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Figure 7-11. Theoretical Variation of Combustion Gain with the 
Stability Parameter N 
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COMPUTES FUEL 
STABILITY DECAY RATE 
CHARACTERISTIC BASED' 
ON O.W.DYKEMA 


I 


SECF25 


COMPUTES OXIDIZER 
STABILITY DECAY RATE 
CHARACTERISTIC BASED 
ON O.W.DYKEMA 


T 


CONTINUE COMBUSTOR 
DESIGN LOOP 


Figure 7-12. Flow Chart for Subroutines SECF24 
and SECF25 Computing the Decay 
Rates F24(Fuel) and F25 (Oxidizer) 
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A definition of the respective frequency modes is given in the following 
Section 7, 2. 3.1 

After calculating the stability numbers, N , for all ten frequency modes, 
the per cent of maximum combustion gain is found from the relationship presente 
in Figure 7-11 which is coded into the program (Subroutine DYKEMA) . 

Based on the assumption that a decay rate of 100 is assigned to be 
equivalent to a value of 0 for the per cent of maximum combustion gain, and 
a growth rate of 100 is equivalent to 100 per cent of maximum gain, then 

F24, F25 = 2 (max. gain - 50) (7-18) 

If the 'fuel or oxidizer is in a gaseous state (ISTATF or ISTATX = 1), Subroutine 
SECF24 or SECF25 will be bypassed, • , . 

Figure 7-12 depicts the flow chart' for the controlling subroutines SECF24 
and SECF25 . , 

7,2.3. 1 Definition of Frequency Modes for Subroutines SECF24 and SECF25 
Without Baffles Present. (Computed in Subroutines ACUSTC) , Figure 7-13 
defines the symbols used in the following equations. 



LI = chamber length without converging section 
LC = chamber length including " " 

LDS= acoustical length of chamber 
D1 = chamber diameter (near injector face) 


Figure 7-13. Combustion Chamber Geometry Definition 
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• LONGITUDINAL MODES . Longitudinal vibrations can be "caused by either ' 

1. A resonance effect analogous to a stretched string or an organ 
pipe where a supposedly homogeneous gas column exercises acoustical 
vibrations of various modes, or by, 

2, Changing gas properties in the combustion chamber (e.g., shifting 
mixtures ratio) generating entropy waves. 

To keep the program input requirements to a minimum only the conventional 
acoustical vibrations (Case 1) are being consi'dered. . ^ 

It can be supposed that the near throat area of the nozzle acts as a 
closed boundary to the pressure f luctuatiojis , Reference . 7-6 . Then a pressure 
pulsation generated at the injector face might propagate downstream' at . the 
speed of sound, be reflected near the throat area, and then return .to_ th,e . 
injector at the speed of sound. The frequency is determined by the time ' 
required for the pressure pulse to .perform the complete cycle. Thus for the 
1st mode 

C 

[freq. 1]: . = , cps (7-19). 

^ ^ 2*LDS 

where C is the speed of sound of the gas in the combustion chamber and 
is defined by 


C =ygkRT . (7-20) 

where 

g = gravitational constant _ 

k = ratio of specific heat 

■ R = gas constant 

T = combustion temperature ... 

When using th.e analogy with the vibrating string, i.e. , that a nodal point at 
the reflection areas exists, the second arid third modes can be represented by 
the following sketch. 
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Accordingly., the second mode of. the longitudinal vibration can be expressed by: 

C - ■ ' ■ 

, cps, , .(7-21) 


[freq. 2] = 


IDS 


and the third mode 


[freq.- 3] : f 


L(3) 


2*LDS 


•; cps* 


(7-2-2) 


Using .the material presented in References 7-7 and 7-8 and considering the 
above mentioned requirements , .tangential^ radial , and combined modes can-be 
presented as follows: 


TANGENTIAL, MODES 


First mode : ' 


[freq. 4] f. 


T(l) 


0.586——, cps (7-23) 
D1 


Second mode: 


[freq. 5]: f^ 


( 2 ) ,• 


0..972 , cps (7-24), 

' ■ ■ D1 ■ ' - 


Third mode: 


[freq. 6] : f. 


T(3) 


1.337 , cps (7-25) 

D1 


RADIAL MODES 


First mode: 


[freq. 7]: .= 1.22 —, cps 


D1 


(7-26) 


Second mode: 


[freq. 8]': f ' ^ = 2.233 —, cps (7-27)' 

.• / . ■ ■■ • ■ Dl 


LONGITUDINAL MODE 1 WITH TANGENTIAL MODE 1 


^LCl), T(l) = V'mI) * *T(1) 
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LONGITUDINAL MODE 1 WITH TANGENTIAL MODE 2 


[freq. 10] : 


^L(l) ,T(2) 



+ 

T(2) 



(7-29) 


7. 2. 3. 2 Definition of Frequency Modes for Subroutines SECF24 and SECF25 
with Baffles Present . If there are radial and/or circumferential baffles 
present, the equation for the frequency, f, that gives maximum value for the 
combustion gain with respect to transverse and radial modes , becomes 


where 


^*^vh 

r*2iT 


(7-30) 


a ■= speed of sound in combustion gas, in/sec 

r = combustion chamber radius near the injector 
plane, ft 

s^h ~ Eigenvalues given in Table 7-7 (Section 7.8.3) 
for various v and h numbers 

V = number of transverse modes 

h = number of radial modes 


The program coding selects the original three longitudinal modes as defined 
in Section 7.2.3. 1 (Subroutine ACUSTIC) and, via Subroutine FNDSNH (see Section 
7.8.3), the lowest six svh-values , satisfying the exclusion rules as given 
in Section 7.8.3 (for more details for program provisions to accept baffles 
see Section 7.8.3). 


7.2.4 Subroutine SECF26. Subroutine SECF26 computes. the decay rate for 
the longitudinal modes of high frequency oscillations based on the sensitive 
time lag model of L. Crocco and S. I. Cheng, Reference 7-3, with the correla- 
tion equations for the interaction index, n, and time lag,T , developed by 
D. T. Harrje and F. H. Reardon, Reference 7-9. According to the Equation 
(3.01.20) of Reference 7-3, the decay rate, F26, can be calculated by *) 

*)This a shorter version of Equation (3.01.20) of Reference 7-3 in accordance 
with the "Special Instructions for Task II" of Contract NAS 3-13331. 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF24 AND SECF 25 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF24 AND SECF25 
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(7-30) 


1 - B exp[2s] 

: — — - - M[(l-vn) + TnexpC-sx)] = F(s) =0 

1 + B exp [2s] . 

where . , ^ ^ 

n 1 1/2Cy-1)M 

“ “ 1 - 1/2(y-1)M 

M. = Mach number of gases in chamber 

s = Lap lace, transformation variable 

_ r 

T = T * 

, 1000 * LDUM 

T = Time lag 

Cg = Speed of sound in' chamber 
LDUM = Acoustical length of chamber 
Y = Ratio of specific heat 
n = Interaction index 

For non-hypergolic propellants with coaxial injection: 

Tj'= 0.076* Fp ,millisec 


(7-31) 


(7-32) 


where 


P Tc/'’crlt " ' ' Cc/Pcrit) 

’’ Vc/Pcrit >1 = 1 

= 1.446 - 0.4*VR*sin(|)+0.0685*(VR*sin4))^ 

yj^ _ Velocity of outer ring propellant 
Velocity of inner hole propellant 

(j) = Impingement aijgle 

n = 0.45 (can be modified by input in COMMIS 
where it is designated NTAUNl) 
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NOTE: T-j is computed in Subroutine TAUNl.. 

VR is computed if input for IC0AX=1, then VR=VFUEL/VOXID; 
‘ . or for IC0AX=2, then VR=VOXID/VFUEL 

For non-hypergolic propellants with non-coaxial injection: 

T2= millisec (7-33) 

where 

= injector orifice diameter 

The interaction index, n, satisfies a fourth order polynomial curve- fit, 
and is 

n = exp[AQ+A^JlnD^ + A2(£nD^)^ + A^C^nDj^)^ + A^CAilD^)^] 

with 


Aq = -55.5352E-02 
Aj = -96.4330E-03 ^ 

A2 = -28.1489E-03 
A3 = -59.2126E-05 
A4 = 66.5459E-04 

These data are coded in Subroutine TAUN2, and .T2 is computed in the same 
subroutine. 

For hypergolic propellants: 

T3 = millisec (7-34) 

where 


= chamber pressure 

n = .7 (can be modified in Subroutine COMMIS and 
is designated NTAUN3) 

T3 is computed in Subroutine TAUN 3, where equals DFVEL or DXVEL, 

Figure 7-14 depicts the structure of the controlling Subroutine SECF26. 

The optimization program AESOP (Subroutine CTLOPT) is used in an inner loop as 
a root finder by minimization of |F(s)jin Equation (7-30). The parameters are 
the elements of the complex root defining decay rate and frequency of oscillation. 
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COMBUSTOR DESIGN 
LOOP 


CTLOPT 
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PROGRAM 
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Figure 7- 14, Flow Chart for Subroutine SECF26 Computing the 
Decay Rate Characteristic for Longitudinal 
Mode, F26 
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V) 




(D 



0 




o 

O 



p 



•H 

4^ 

i-H 






T3 

U 



0 

O 



4-> 

4-> 

•H +J 


3 

O 

'O V) 


P. 

0 



g 

*r-i 

+-> f-f 


O 

P 

: «H O 


U 

•H 

‘ W) 4-> 
C V) 



X 

o 


(/) 

P 

- 1-t x> 


o3 


£ 

c 

bO 

P 

^ o 

o 


0 

O o 

•H 

P 

X 

4-1 

4-> 

O 

g 

CD 

O 

•H 

p 

4-> bO 

0 


X 

C 03 


C/) 

u 

03 ?-( C 

C ' 

P 


+J CD O 

•H 

XI 

p 

C/) > -H 


g 

o 

C 03 +J 

4^ 

o 

• H 

O o3 

O 

u 

X 

O CD 



C/) 

x: a- 

0 

p ^ 

P 

^ -M <D 

i-H 

•H 

X 

O 

bO 

H 

g 

4-> C O 

C 

'O CO 

o 

O O *H 

o3 

P o 

o 

03 +-> 


p p. 

44 

4-1 »-* C/) 

4-> 

o 

o3 *H 

c 

t/1 0 

o 

bO «H H 

0 

P 


C 4-> CD 

g 

4h *H 

p 

•H 4-* 

0 

O +-» 

0 

4-> CD O 

bO 

P 

X 

x: C 03 

C . 

O 

0 

bO O 

•H 

0 U 

g 

•H pt, o3 

& 

0 XI 

P 

CD X x: 

g 

p. p 

•H 

^ CD CJ 

»— 1 

c/) cn 

Q 




b 0 CI 4 

hH 

1— 1 

0 

I-H 

• ^ 

HH 0 





X 


0 

0 

0 

< UJ 

r-H 

00 

> P 

> p 

Z X 


1—4 

•H 0 

•H 0 

1-4 X 

P 


X X 

X X 

S 

0 

0 

-H ^ 

•H g 

0 > 


X 

V) n 

V) p 

z 

0 


0 :z: 

0 :z: 



0 

CL 

CL 



0 

0 

t/> 

CO 


0 

0 

0 

H 

I-H 

1 



-g 

§ 


0 

P^ 

P 


T3 

•H 

• H 

X 

< X 
CD 0 

H-4 OQ 



CO 

u 

0 

Q 

3 ^ 
b > 

j 

-e- 

CJ CO 
tij 










^ X 

vD 

Cl 



0 

OJ 

S 



0 OQ 

(X 

h - 1 

CO 

CD 

0 S 
a: > 

CL CO 

< 

0 

0 

Q 



DFVEL D. inches Positive I Injector orifice diameter for fuel or oxidizer. For 

or ^ Number definition see Section 2.4 

DXVEL 


(/) 

0> 

(/) 

(/) 

(D 


4-» 


CO 


0 

0 


0 

TD 


P 

0 



g 


0 



X 

I-H 


X 

P 



P 


X 

•H 


0 

T3 



P 

X 

X 

X 

p 

p 

■H 

p 

p 

W) 

X 

p. 

P 

CO 


0 

p. 

X 

i-H 

0 

p 


0 

p 

P 


p 

0 

I-H 

•H 

X 

p 

bO 


p 

P 

0 • 

0 

g 

X 

•H 


p 

X 


p 

p 

P 


X 

0 

X 

•H 


p 

> 

I-H 

0 

P 

P 

0 

P 

0 

Q 

CD 

CC 


0 

1 




0 

1 



> 

P 


p 

CN 

> 

P 



•H 

bJO 


0 

I-H 

•H 

bO 


0 

X 

0 




X 

0 


P 

•H 

12: 

0 

g 

CNJ 

•H 

:z; 

0 

rX 

to 


> 

P 

• 

to 


> 

P 

0 

P 

•rH 


CM 

0 

p 

•H 

> 

CL 

0 

X 


tr> 

Cl, 

0 

X 



O 

<D 

I 

c 

•H 


] 


too 




cti 


I-H 


PL 

vO 


VO 

\o 

W 

CM 

P 

CM 

(N 

cc 

PL 

0 

PL. 

PU 

CP 

Di 


a: 



CP 


CP 


205 


ICOAX — — 0, 1, or I CoaxiAl injection indicator. For definition see 

2 Section 2.4 












PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF26 
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.Mach number in combustion chamber (computed in Sub- 
routine POSTAL) 












PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF26 
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PARANETERS (COMMON) USED IN SUBROUTINE (S) SECF26 
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When a root is found, it is swept out by a root sweeping routine, SWEEPR,; 
and another root' solution is attempted. This process is continued for a 
predetermined number of roots. The resulting complex roots are examined, and, . 
from a series of solutions possible, the roots with the largest positive real 
part having the three smallest values of the non -negative imaginary part 
are selected. 

The characteristic function F26 is the real part of the root multiplied 
by a factor which considers the combustion chamber geometry. It becomes 

CS 

F26 = Re * (7-35) 

LDS 

CS = speed of sound in combustion chamber 

f 

LDS = acoustical length of combustion chamber 


7.2.5 Subroutine SECF27 . Subroutine SECF2 7 computes the stability 
decay rate for the transverse modes of high frequency oscillations, F27, 
based on the sensitive time lag model of L. Crocco, D. T. Harrje, and F. H. 
Reardon, Reference 7-7, with the correlation equations for the interaction 
index, n, and time lag, t, developed by D. T. Harrje and F. H. Reardon, 
Reference 7-9. The interaction index, n, is twice that of Subroutine SECF26, 
and the time lag T...,is the same as. in Subroutine SECF26. SECF27 considers the 
boundary condition at the convergent nozzle entrance (nozzle admittance) in an. 
approximate manner by the complex admittance coefficient E (frequency-, mode and 
nozzle entrance velocity) , as well as by the geometric relationship of the 
combustion chamber, RA. 

The subroutine for solution of F27 is similar in structure to that 
employed for the computation of F26. The same routines TAUNl, TAUN2, and 
TAUN3 are utilized in the computation. 

Figure 7-15 shows the structure of the controlling Sibroutine SECF27. 

The primary difference between this routine and the controlling subroutine 
SECF26 is in the stability characteristic equation employed. This equation is 
coded in Subroutine TCOMBI. SECF27 calls subroutine TCOMBI on an inner AESOP 
loop. TCOMBI generates the complex function for neutral stability 

hj P + h2 = F(s) SO (7-36) 

as defined in Equation (28) of Reference (7-7). In this equation is: 


where 
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TAUfJl 


SENSITIVE TIME LAG 
FOR NON-HYPERGOLIC 
PROPELLANTS 

(coaxial injection) 


TAUH2 


SENSITIVE TIME LAG 
FOR NON-HYPERGOLIC 
PROPELLANTS 

(non-coaxial 

INJECTION) 


♦ 


TAUN3 


SENSITIVE TIME LAG 
FOR HYPERGOLIC 
PROPELLANTS 


BDATA7 


INITIALIZES ALL DATA 
IN INTERNAL AESOP 
DATA BASE AND TURNS 
OFF OUTPUT INDICATORS 

4 — 


COKibUSTOR DESIGN 
LOOP 



SUBROUTINE TO CALL 
IN THE RESPECTIVE 
SENSITIVE TIME LAG 
taunI, TAUN2, OR 
taun3 


CTLOPT 


OPTIMIZER CONTROL 
SUBROUTINE FOR AESOP 
PROGRAM 


TCOnBI 


COMPLEX CHARACTERIS- 
TIC EQUATION FOR 
NEUTRAL OSCILLATIONS 
AT THE STABILITY 
LIMITS 


SWEEPR 


ELIMINATES PRE- 
I VIOUSLY LOCATED 
ROOTS 


bNDCHK 


DETECTS BOUNDARY 
VALUE CONDITIONS 
(root OUTSIDE OF' 
REGION OF interest) 


i 


SECF27 


COMPUTES DECAY RATE 
CHARACTERISTICS FOR 
TRANSVERSE MODES 
BASED ON CROCCO AND 
CHENG 


T 


CONTINUE COMBUSTOR 
DESIGN LOOP 


FNDSNH 


DETERMINES THE 
RESPECTIVE Sj;/,- 
VALUES FOR A GIVEN 
BAFFLE DESIGN 


Figure 7-15. Flow Chart for Subroutine SECF27 Computing the 

Sensitive Time Lag Characteristic for Transverse 
Modes, F27 
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and 




hi = ■YUe[l-is^l^EJ (u/Ug)-dz] 
0 

P = n(l-e“^^) 


(7-37) 


, (7-38) 


and, omitting parts pf the function in accordance with the Spepial Instructions 
for Task II of Contract NAS 3-I333I, 

.1 . _ 

. 'h2 = - i(f --^svhze 

- (f--) + i(Y+l) s^u f (^)dz] 

f J 2 f ^ J ^e 

' 0 (7-39) 


and 


I-Y to u 


E = - 


(7-40) 


2s, c’ 
vh 0 


Note; 


The expression for E is derived from Equation (26a) in 
Reference 7-7 as follows: 


E = 


-fP yU^ 
0 ' e 


P P 
e o 


= -Yf— 


(7-41) 


Introducing f = to/s^h 


U 


E =.^ 


^;h p 


(7-42) 


Furthermore, assume that the nozzle is distributed nozzle 
so that the exit velocity of the chamber is at a constant 
Mach number then 

U* = — ^yRT* 
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CT*/Tq) u/Co 


Furthermore, assuming adiabatic law 


Y-1 



Therefore • 

U = P (7-43) 

e e Co 

Using a linear analysis for Equation (7-43) 


^e 



(7-44) 


Substituting Equation (7-44) in Equation (7-42) 


-yo) y-1 u 
^h 2y Cq 


= ill 
2 


to 


’vh 



(7-45) 
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Symbols in Equations (7-36) through (7-45) mean: 

Cq = Speed of sound in combustion chamber 

E = Nozzle admittance coefficient 

f = Reduced frequency = 

n = Interaction index of propellants 

P = Proportionality factor relating burning rate 
perturbation to pressure perturbation 

P^ = Combustion chamber pressure at nozzle entrance 

P = Combustion chamber pressure near injector 
o 

P = Dimensional quantity of the chamber pressure 

R = Gas constant of the combustion gas 

s = Laplace transformation variable 

s = Mode of oscillation constant 
vh 

= Combustion temperature .near injector 

T* = Dimensional quantity of combustion temperature 

U = Velocity perturbation distribution function at 
nozzle entrance conditions 

"k 

U = Dimensional quantity of the velocity perturbation 
distribution function 

u = Steady state axial gas velocity 

Ug = Steady state chamber exit velocity 

2 g = Chamber length to exit 

Y = Specific heat ratio 

0 ) = Angular frequency 

T = Time lag 
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In accordance with NASA's technical monitor, the expressions for h^ in 
Equation (7-37) and h 2 in Equation (7-39) are modified for coding purposes as 
follows 

(1) The integral expression 
geometric relationship of the combustion chamber 

RA = LC/R 

= 4LC/(DI + D2) (7-46) 

(2) Zg becomes RA 

(3) Ug becomes Mach number M ' 
where 

LC = Combustion chamber length 

R = Mean radius of combustion chamber 

D1 = Combustion chamber diameter at injector face 

D2 = Combustion chamber diameter at the onset of 
the converging section 

Thus 

(7-47) 


(7-48) 


where M is the mean Mach number in the combustion chamber. 

The optimization program AESOP is used as a root finder by minimization 
of JF(sJ I in Equation (7-36) for the various modes of oscillations. The 
optimization parameters are the elements of the position in the complex plane 
The resulting roots define decay rate and frequency of oscillation. When 
a root is found it is swept out by the root sweeping routine, SWEEPR, and 
another root solution is attempted. This process is continued for a pre- 
determined number of roots. The characteristic function F2 7 is the real 
part of this root multiplied by a factor which considers the gebraetry of the 
combustion chamber: 


hj = YM(l-i-s^j^ E * RA) , . 

h2 = -(Y+l)M-i(f-^s^j^ RA + 
ra2 


Ze 

is being replaced by the 
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where 


• F27 = Re(s) * 


CS = speed of sound in combustion chamber 
DC = Combustion chamber diameter 

A characteristic function for both, fuel and oxidizer, is determined in this 
manner. Currently, the larger of the two yalues is retained. 

At the present time the critical value of F27 is obtained for the seven 
modes listed in Reference 7-7, and presented in Table 7-3. Within the program, 
the frequency modes are stored in Subroutine ,CMBINI as SNUH(l) through SNUH(7) . 



Purely Tangential 
Modes 

Purely Radial 
Modes 

Combined 

Modes 

First 

Second 

First 

Second 

- 

- 

- 

II 

1.84129 

3.0543 

3.8317 

7.0156 ■■ 

5.3313 

8.52'63 

6.7060 


TABLE 7-3. Modes of Oscillation for Transverse 
Combustion Instability in a’ 

Circular Chamber 

t 




7. 2. 5.1 Definition of Frequency Modes for Subroutine SECF27 with Baffles 
Present . If th^re are radial’ and/or circumferential baffles present,, the program 
selects Svh'values determined by Subroutine FNDSNH under' the condition’ that s^j^- 
values have to be omitted according to; 
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V ^ F81/2, if F81 is an even number 

V i F81, if F81 is an odd number 


and 


h > F82, if F82 is an odd or an even number. 
Furthermore, two additional conditions have to be met: 


( 1 ) 


The integer number, Uq, of values has to satisfy 

18 




s^h(min) 


C7-50) 


where Si^hrminl lowest allowable s^j^-value selected from 

Table 7-4, Section 7.8.3 


C2) 


s^j^ values have to satisfy 


vh 


r*6280 

a*T 


(7-51) 


where 

r = combustion chamber radius at the injector plane, ft 

a = speed of sound in combustion gas, ft/sec 

T = sensitive time lag as used before in Subroutine 
SECF26, millisec 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF27 
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0<MC<1 0 Mach number in combustion chamber (computed in 

Subroutine POSTAL) ■ 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SE.CF27 
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7.2.6 Subroutine SECF28 . Subroutine SECF28 computes the stability 
decay rate for various modes of high-frequency oscillations (longitudinal, 
transverse, radial), F28; it is based on the acoustic wave solutions for 
combustion response by R.J.Priem and E.J.Rice, Reference 7-10. The 
method uses the flow response functions as determined by Equatibn(ll) in 
Reference 7-11 for liquid propellants, and by Equation(18) in Reference 
7-12 for gaseous propellants. Because system stability is obtained when 
the combustion response line intersects the flow response line at the 
same frequency in the complex plane: 


(Frequency .of combustion response,N ) = (Frequency of flow response, N ) 

C FL 


or 


Nc - = F(s) H 0 


(7-52) 


Specific equations which have been coded for computing the combustion 
response, N^, are Equations (6) , (13) ,■ and (14) of Reference 7-10: 


= Re(N^) + iIm(N(,) 


1 

Y 




-liB. 


1 - 


(B^e + B Ce R ^) 






M[w(e R ^ + Ce + M(B^e ^ + B^Ce ^ ^)] 

(7-53) 


1 - -i-B. 

^ " Cex R 

2 


where 


®1 = 




(M"^-l) (m^-w^) 


1 - M" 


(7-54) 


B2 =. 


t-f 


uM- \/oj^M^+(M''-l) (m^-oj^) 


1 - M" 


(7-55) 


' Bj(1-M^+YGM^) +u(yG-1)M ■ 
. B2(1-M^+YGM^) +w(yG-1)M . 


(7-56) 
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where * . ‘ 

Nj, =■ Combustion response . ■ , ■ 

L = Length of combustion chamber . :• ...■ ' ' ■ 

w ■ R = Radius of combustion chamber at injector face 

M = Mach number in combustion chamber , 

‘ .. '• Y' = Ratio of specific heat of combustion gas 

G = Nozzle response factor 

m = Argument of Bessel function satifying wall boundary 
conditions : 

- > . ■ ' ' ■ The' following values are. used: . 

(1) First transverse mode = 1..8412 

(2) Second transverse mode = 3.0543 

(3) Third trainsverse mode = 4.2012 

(4) First radial mode = 3.8317 

(5) Longitudinal mode =0. 

(The m-values are stored as BESARG(l) to (5) 
in Sabroutine CMBINI) . 

( 1 ) = Angular frequency of oscillations 


To compute the flow response, NpL, it is 


FL 


OF 


OF+1 


■N + 
OX 


— N. 


OF+1 


(7-57) 


where 

OF = Oxidizer to fuel mixture ratio 
Nq^ = Flow response for oxidizer 
Np = Flow response for fuel 
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If the propellant is in a liquid state, the transfer function, Equation 
(11) of Reference 7-11, is used: 


I 

W 

liquid p 

c 





1+2t s 

V 


c T 

l+Cl-23) T s 

A B b 


c T 

1 + -£_JLt ’s 
A B b 


(7-58) 


If the propellant is in a gaseous state, the transfer function. Equation 
(18) of Reference 7-12, is used (modified shorter version according to 
the Special Instructions of Contract NAS 3-13331). Based on the definitions 
depicted in Figure 7-16, it is: 






Aj = Cross-sectional area of injector (injection area) 
b = Vapor^ pressure/ liquid temperature coefficient 

- - J ' , 

Cp = Specific heat of liquid 
g = Gravitational constant 
L = Hydrogen annulus length 
= Combustion. chamber pressure 
P, = Dome pressure '= P- + AP • 

a ■ ■ '' ■ .. 

AP = Injection pressure drop 

Tj^ = Temperature of propellant 

V. = Injector dome volume r 

W =.Flow rate through injector 2 

3 = Vapor pressure-comubstion chamber pressure 
parameter - i ■ 

Y = Specific hbat ratio for c6M)ustibn'’ga§' ' 

YH = Specific heat, ratio for hydrogen 

^ X = Latent heat of vaporization 

; p = Hydrogen density , 

T|j = Time delay constant for burning 

T|jS = Phase angle due to delay in burning =1.9 rad 

s = Laplace operator ' 

Ty = Mean drop. life, time = M/W (coded as TAUV = 
L50/VINJ*2*CS/DC in Subroutine FRPROP) 


Figure 7-17 depicts the structure of the controlling Subroutine SECF28. 

The equation for combustion response is computed in Subroutine ACWAVE, and the 
equation.. for >fclbw response is computed in Subroutine FLORES. This subroutine 
also calls upon subroutine FRPROP to calculate the flow response for individual 
propellants (using Equation (11) of Reference 7-11, and Equation (18) 
of Reference 7-12, respectively). 
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Figure 7-17. Flow Chart for Subroutine SECF28 Computing the 
Decay Rate for Various Modes of High-Frequency 
- . Oscillations, F28 • 
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The optimization program AESOP is here also used as a root finder by 
minimizing of |F(s)| in Equation (7-50). The optimization parameters are 
the elements of position in the complex plane. The resulting complex roots 
define decay rate and frequency of oscillation. The frequencies of both 
response functions are compared in Subroutine CTLOPT which finds the solution 
by driving differences to zero for the various Bessel arguments, m. 

The resulting characteristic function F28 is then the real part of that 
root which is least stable, multiplied by a factor which depends on the 
geometry of the combustion chamber; 

F28 = Re(s)* (7-62) 

DC/2 


where 


CS = Speed of sound in combustion gas 
DC = Combustion chamber diameter 


7 . 2 . 6 . 1 Definition of Frequency Modes for Subroutine SECF28 with 
Baffles Present^ If there are radial and/or circumferential baffles present , 
the values of the Bessel arguments which appear in the equations for computing 
the combustion response (Equations (7-54) and (7-55)), have been replaced by 
^Vh values, except for the case of m=0 . The program selects the lowest four 

s , values via Subroutine FNDSNH which are not forbidden by 

Vh . .. 

V i F81/2, if F81 is an even number 

V - F81, if F81 is an odd number 


and 


h > F82, if F82 is an odd or an even number. 


Furthermore, the following condition has to be met; the integer number 
of values has to satisfy 

4 ^ = § (7-63) 

^vh (min) 


where s,,i, is the lowest allowable 

Subrout 


s^h value 


selected from Table 7-4 via 


In addition, the longitudinal mode (m=0) of the original code is retained. 


228 



PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF28 










U 











0 




. 





/ — S 


M 




i-H 





. \0 


•H 




• 




LO 

1 






r^ 




rH ^ 

(N 


•H 


0 




0 


.-f 0 



X 


a 


c 

+-> 

0 


Cd cx: 

0 


0 

p 

p 


o 

V) .H 

to 


^ < 




•iH 

44 


•H 

0 r-l 



CO 

P 


?H 




+-> 

T3 .H 



bO ^ 

bo 


0 


p 


o 

o ^ 

. 


q cQ 


f N 


0 

0 


(D 


0 


•H 

Uu 

X 

I-H 

44 

44 


CO 

. +j 

4 -> 


X ^ 


1 

0 

P 

0 



TJ (/) 

cd 


44 bO 

0 

CM 

P 

&• 

0 


0 ) 

0 

?H 


to 0 

0 

r 

44 

£ 

•r-i 


<D 

c ^ 



•H 0 

to 

0 

V / 

0 

P 


in 

•H O 

X 


+4 U 

V — / 

U 


0 

' *H 



+-> 

cd 


Cd X 


P 

44 




• »\ 

e V) 

0 


to +4 

to 

bO 

c 


44 


c 

o :3 

0 



0 

•H 

p 

to 

P 

o 

o 


'O 


C 

I-H 

Pu 

I-H 

P 


hH 

•H 




0 i-H HH 

44 


I-H 

bO 

P 

H 

+-> 

^-1 o 

4-1 


•H 

44 

0 

0 


0 

cu 


o o 

0 


+4 e? i-t 

p 

0 

PL, 

P 

X 

hH 


4 h 



0 Oi OQ 

X 

to 

0 

0 

B 

OC 

•H 

0 

0 


c < s 


V — * 

U 

• H 

P 

u 

U 

4-) bO 

W) 


P CO CJ 

I-H 


Ph 

44 

X 

CO 

+J 

cs 

P 


44 m 

P 

to 


to 

a 

Ol 

w 

^ C 

cd 


CQ 0 

•H 

0 


P 


Q 

• H 

4-> 0 O 



•-H C 

+4 

r— 4 

• iH 

X 

p 



W > -H 



0 -H 

P 

44 

P 

£ 

0 



C oi 4-> 

X 


to • +4 

0 

44 

cr 

0 

•H 



o oj 

u 


to to p 


p 

• H 

0 

44 


<D 

o o ;3 

M 


0 p 0 

0 

X 

I-H 

/■ — \ 

to 


rO 

X o" 

cd 


OQ 0 H 

44 



P i-J 

P 



^1 4-> 0 

0 


•H X 

S 

I-H 

44 

•H < 

X 


3 

O 

to 


0 +4 P 

P 

P 

0 


£ 


(i: 

+-> C u 



X *H CO 

CJ 

•H 


T3 CO 

0 



O O-H 


(N 

44 T3 

JH 

T3 

44 

P 0 

0 


j-f 


0 

• 

P c 

• H 

P 

P 

P o. 



o 

44 I-H to 

44 

LO 

44 0 *H 

a 


0 

0 

44 


4^ 

cd t4 


• 

0 a 



X 

to 0 

0 



bO«H fH. 

4-> 

CN 

T3 

44 

44 


p 



+J 

4-> 0 

P 


44 X 0 

0 

0 

0 

44 -H 

P 


P 

•H C 4-> 

0 

P 

C ^ 



•H 

0 44 

0 


oJ 

+J 0 O 

B 

0 

0 cd 0 

u 

?H 

44 

P 

44 


+-> 

tJC C cd 

0 

•H 

e T) +4 

0 

0 

•H 

TJ 0 

0 


(/) 

W) o H 

4J 

+4 

P P to 

X 


0 

0 Pi 

B 


C 

•H pH cd 

Cd 

0 

bO p 

B 

£ 

0 


P 


O 

^ XX 

+4 

0 

H 0 to 

P 

p 

pH 

PP P 

•H 


u 

^ 013 

CO 

CO 

< X -H 

'Z 

z 

CO 

CO CO 

Q 


:d o a- 

1-1 O 


X 








0 

P 

0 

p 


0 

0 

0 

< 

w 

i“H 


•I 




> 

0 

> 

0 


> 

> 

> 

z 

X 0 



CM 

to 

CM 


•H 

•bO 

•H 

bO 


•H 0 

•H 0 

•H 0 

I-H 

X 

P 

0 

I-H 


t-H 

I-H 


0 0 


0 

0 

44 p 

X P 

X P 

S 

«— 1 

0 

to 


LO 

0 

to 

•H 

44 

•H 

X 


•H »-H 

•H »-H 

•H X 

0 

> 



00 

0 

CM 

00 0 

to 

P P 

to 

P 

P 

to P 

to P 

to P 

. z 


0 


• 

• 

• 

• • 

0 

HH 0 

0 

HH 

0 

0 > 

0 > 

0 > 





I-H 

to 

''d- 

to 0 

Cm 


Pu 



CU 

PU 

(X 














oc 

0 


to 


CO 












B 

0 

0 


H 1 









1 



^ X 

0 

X 


5 ' 

1 

1 

1 




1 


i 



X 


0 

p 


5 












p 

c 

•H 














CQ 

•H 


X 

< 

X 















u 

0 















I-H 

S 












Ph 

0 

u 

§ 

$ ^ 

1 

1 

B 




1 


1 



0 

U 

Q 


^ 0 



00 

C3 


Q 



CO 

00 

CM 


I-H 

< 

Ph 

X 

IS> CQ 

I-H 

CM 

§ 

< 

u 


Oh 

P Oh 

0 s: 

Q 

Uh 

CO 

X 

X 

U 

0 u 

OC X 

< 

< 

X 

w 

Ph 

pH 



Q- CO 



< 

CQ 

CQ 

CQ 




229 








Q> 


o 

<D 1 




1— 1 
Z 

CM 

O 1 




> 


> 

> cd 




HH 

i-H 

> Cd 


• 

• 


<u 

•rt O 

•H W) 


<D 


VO 00 

•k 

•H bo 

0 

m 

w 

4-> 


'P 

P o 


;3 

Oi 

VO S 

<N 

4_) 0) 

:3 

CO 

p 

•H 


•H 

•H :z: 

0> 



»-• u 

• 

•H Z 

rH 

(-0 

O :Oa. 

V) 

cd 

(/) cd 

(/) 

> 

cd 

t—i 

Ol 

CM 

CO > 

Cd 

< 


o 

> 

O > 

o u 

•H 

> 


• C 

to 

O P T-l 

> 

pu 

• 

Oi 


Ci. 

Oh O 

P 



•H 


Cu O P 





to 

0 

CNJ 


(N 


jC 

10 

• 

o 


: a 

4-T^ 

T3 

Cd 

0 


•H 

Xi 

i— 1 

P 

•H 




t-J pH 

W 5-1 W U Cu 00 

> O > hJ O ^ CVJ 

u. X W W P-. 

Q Q a Q 


CO 

00 

ou 

CO 


oc 

00 


HH 

E-< 

w * 

pH 

CM 


00 

<x> 

Cti * 

PJ 

PU 

p 

CM 

CM 

z 

02 

PJ 

o 

pH 

PL, 

C3 

c:? 

cc: 


pc: 




o 











PARAMETERS (COMMON); USED IN SUBROUTINE (S) SECF28 






X 

nd 



t 



} 




:z 











CO 

13 










Q 

" lo or* 




P 







03 -H ' 



U' 

0 






PU 

bO 

0 


0 

P 








. 4 :: 


M 

o 


CD 




0 

II II 

P 


•H 

0 


CD 




P 




'O' 

•r^ 


C/) 




•H 

« (N 

‘ P 


•H 

c 



CO 

O- 

00 

P 


•H 


X 

•H 


p • 

PU 

o 

(N 

P 




0 



o 

C3^ 

aC 

PP 

O 


nd 



B 



o 

Cu 

o 

P 

• P 

0 


p 

0 




cc: 

w 


•H 



0 

P 

■; 

•H 


pu- 

CO 

■ p 

X 

0 



Mu' 


P 




CO 

• -0 

0 


fp 



■H 

(D 

CD . 

CD. 





0 

0 



P 

P 

P 

s 

p 

to 


p 

0 


(D 

•H 

• H 

•H 

o 

o 

p 


M-l 




P 

P - 

P, 

p 

p 

PI 



..TO 



:3 

P ' 

P 

4h 

' o 

0 


P 

^.P 


P 

o 

O 

O 


0 

S 


0 

. to 

> 

O 

u 

u 

P 

to 

•r-» 



4h ‘ 

• H 


4 -» 




0 


bO 



'Id 



p 


P 

P 

■ »H 

■ P 


IP 

'■ — •' 


• 

CO 


CO:; 

Pt 


CtJ 


0 



X 




oi 

P 



' "'H ' 

4P 


CD 


P 

P 

> 

ctS 

r— < 


P 

P 

* *■ - i 

T? 

o 

o 

O 



0 


' o3 

bO 


C 

tp 


4 -i 

rC 

0 

to 


tsi 

P 


•H 





P 

to 


U »H 

0 

: 0 


c 

C 

c. 

to 

cS 

0 


^ P 

iP ^ N 

•H 

X 

o 

O 

o 


P 

CQ 


■ 0 

TO 


p 

•iH 

•H 

•H 


' to 



^ Ph 

p 0 


•H 

p 

P 

p 

0 to 


tP 


TO 

0 p 

. P " . 

X 

Ph 

P 4 

a 

p • 

p 

. 0 



X TO 

0 * 


o 

O 

o 

O 00 

■ 0 



■ i ■< 

6 

’ , 

. 

f— < 




0 * 


P 


tp 

TO P 

' iP- 

Ph 


CD 

CD 


•H 

0 


• 0 

X TO 

0 ' 

6 

p 

P 

P 

0 

nd 



t j ■ 

0 0 

. P '■ •■ 

O ^ 

p 

•H 

•H 

to p 

•P 

‘ B 


' P 

P 

. ; 

u . 

o 

P 

P 

o 

X 

P 


TO' 

P 4P 


CM 


• 


Cp -H ^ 

o 

■ C 


0 

0 P 


u ' 

bfi 



O P 




-■ rP 

•H 

■ ' 

o p 


P 

P 

CJ 

' p 

g 

s 


P 0 


P o 

•H 

O 

O 

. p 0 

* o 

. P 

TO 

P 

to P 

4 :: ■ ^ 

O «H 

p 



0 CO 


B 

P 

p: 

P 

p. 

CD P. 

•H 

nd ' . 

nd 

X 

i-H 

•H 

bo 

' . 0 

4D 0 

bio 

•t-^ O 

P 

oS . 

ccJ 

B o 

0 

/ ^ 

0 

' p 

B 0 

P" 

c o 


CD 

0 . 

p 0 

. P : 

'• 03 

P 

’ TO 

0 TO 

0 

l-H CO 


ct: 

Qc:. 

Z to 

*. Uh 

• S 

PP 

- nJ 

{ 

U mh 

pj - 



p 


p‘ 

P 








.0 


p 

0 




0 0 

CD 

0 J 







CM 


> 

» ^ 

>.- 

9S 

»\ 


•w 


0 


. 

•H 0 

•H 0 ^ 

•H 0? 

fH 

TtJ- to 

tp 

IP (N 

rP CM 

P 

P 

0 

p TO 

P P 

P TOt 







0 

P 

•H P 0 

•HP 

•H P; 


•S 


. 

•> 




' to TO 

to TO ■ 

to TO 

. 0 

to 

0 

0 

0 


P 


- 0 > 

b > 

o.>- 









cu 

cu 

CU 









B 

to 

to 









4D ■ 

0 

0 

1 


1 

I 

1’ 

! 

I 

: I 


•s : 

. -s . 









P 

•' P • 

p 









P 

- »H 

•H 









CQ 




231 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF28 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF28 
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Positive 0 Fuel or oxidizer mass flow rate through injector; 

Value see Section 7.9 . 















7.2.7 Subroutine SECF29 . Subroutine SECF29 computes the decay rate, 
F29, based on the non-linear stability analysis of R. J. Priem and D. C. 
Guentert, Reference 7-13. This analysis calculates the regions of high- 
.frequency combustion instability from two models which determine the local 
burning rate. The stability decay rate becomes 

F29 = -1001ogiQ(Ap/.04) (7-64) 

where Ap is the pressure amplitude and is obtained from Figure 7-18 as a 
function of the buming-rate parameter, L, and the reduced velocity 
'difference AV. 

In Figure 7-18 

L = Buming-rate parameter 


where 

= Diameter of combustion 

^ “ Diameter of combustion 

injector plane 

M = . 25 /L 50 

■ LgQ = Length of combustion chamber required to 
vaporize 50 per cent of the propellant as • 
determined in subroutines SECF20. and SECF21 

= Considers the presence of radial baffles as 
follows*: 

• ' (1) becomes a positive 'integer and is 

equal to F81, if F81 is an odd integer, and . 

(2) becomes a positive integer and is equal , 

to F81/2, if F81 is an even integer ' 

If no baffles are present, Myy becomes 1, and 
thus the original expression for the burning-rate 
parameter, L, is retained. The expression for L 
is coded is subroutine NLSTAB, for which subroutine 
SECF29 is the controlling subroutine 


♦According to Reference 7-13, Appendix B, the radial velocity and all deriva- 
tives in the radial direction are assumed to be zero and, therefore, the effect 
of circumferential baffles, is omitted. 


M 


°C ^XX 


chamber throat 
chamber near the 
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Figure 7-18. Stability Limits of Vaporization-Rate 
and Chemical-Reaction-Rate Models . 
(Figure 8a of Reference 7-13). 






AV = Reduced velocity difference between gases and 
liquid. It is obtained by assuming the gas 
velocity to vary linearly with chamber length, 
liquid velocity is constant, and an average 
velocity is then calculated using a minimum of 
50 fps difference in averaging. It becomes: 


1 (V - 50)V 2 (V„- V,- 50)^v 

AV = — [ 50 + ^ + — i §-] 

a 2(Vj^+ 50 )(Vl+ Vg) 2(Vg'- 


if V^<50, then 

if V -Vt< 50 then 

g L 

where 


Vj^- 50=0; 

V - V - 50=0 
g L. 


(7-65) 


a = Speed of sound in combustion gas 
Vg = Final gas velocity in combustion chamber 
= Liquid velocity of propellant with largest 


The data of Figure 7-18 is stored in Subroutine NLSTAB and_ is treated 
as a three-dimensional surface. The pressure amplitude A (z) is determined 
for a -given burning rate parameter (x) and reduced velocity (y) by the use 
of a general purpose two-dimensional table look-up routine. This subroutine 
incorporates arbitrary scale transformations. In the present application 
log-log scale transformations in both the x-y and y-z planes are employed. 

The burning rate parameter, L, and the reduced velocity difference, AV, 
are also calculated in Subroutine NLSTAB. The coding is done in such a way 
that only the largest Ljq is considered. If the propellant state is a gas 
plus a liquid, then only the liquid is considered. If both propellant 
components are gaseous, then the subroutine is not used at all. 

Figure 7-19 depicts the structure of the controlling Subroutine SECF29. 
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LFIFTY 

COMPUTATION OF 
LENGTH OF VAPOR- 
IZED MASS (detail 
IN FIG. 7-10) 


COMBUSTOR DESIGN 
LOOP 


TLU2 


NLSTAB 1 

, 

■ SECF29 

TABLE LOOK-UP 
ROUTINE 


EXECUTES calculation; 
INCORPORATES THE | 


COMPUTES THE DECAY 
RATE, f29, based ON' 


. . p 

BAFFLE EXCLUSION I 
RULE. AND. STORES DATA] 
FOR TABLE LOOK-UP I 


THE NON-LINEAR 
STABILITY, ANALYSIS 

OF TN d-1409 


CONTINUE COMBUSTOR 
-DESIGN LOOP 


Figure 7-19. Flow Chart for Subroutine SECF29 for 
Calculating the Decay Rate Based on 
Non-Linear Stability Analysis 
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7.3 Pressure Drop Characteristic Modules (Subroutines) 


The Subroutines SECF31 and SECF32 for the injector pressure drop 
'■ characteristics are only output subroutines. The actual computation is 
■ executed in Subroutine POSTAL independently of the weighting factors AF31 
and AF32 which may be equal to zero. However, the injector pressure drop-, 
AP, is required for: 

1. Computation of the chugging stability characteristics F20;.and F21, 
and there for- the following parameters: 

P (pressure drop parameter) 

E (elasticity parameter) 

J (inertia parameter) 

2. Computing the acoustic stability characteristic, F28, and there for 

' the transfer function (Equation (7-59)). 

3. Computing the specific pressure drop characteristics F31 and F32. 


7.3.1 Subroutine SECF31 . Subroutine SECF31 is activated as long as the 
weighting factor AF31 is not zero. The pressure drop for thie fuel side^of 
the injector becomes * 


A P^ — 


2*g*CD^^*Pf(^DF^ NF)2 


(7-66) 


where 


C^ = Discharge coefficient for fuel orifices. 
Df 


DF = Equivalent diameter of fuel hole (smallest- 
diameter) - 

g = Gravitational constant 

NF = Number of fuel holes 


W£ = WT/(1+0F); fuel flow rate through injector 
WT = Total propellant flow rate through injector 
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P£ = Fuel density at injector entry ■ 

OF = Propellant mixture ratio 

Because the computer program is able to take into account a pump-fed engine 
with preburnCTs , the equation for the pressure drop takes the following form: 


WINJF^ 


^Pf = 


IT 


2*g*Cp^^ * pf (jDF‘'*NF) 


(7-67) 


where 

WINJF = WF + WXPREB - WFCOOL 

(For definitions see Figures 2-10 and 2-11 in Section 2.4).- . Finally, 


F31 = Ap. 


(7-68) 


7.3.2 Subroutine SECF32 . Subroutine SECF32 is activated as long as the 
weighting factor AF32 is not zero. The pressure drop for the oxidizer side of 
the injector becomes 


AP 

ox 


iv 2 

03C 

2.g.Cp^".p„^CjDX".N0X)^ 


(7-69) 


where 


C = Discharge coefficient for fuel orifices 
X 

DX = Equivalent diameter of oxidizer hole (smallest 
diameter) 

g = Gravitational constant 
NOX = Number of oxidizer holes 

W^x WT/(1+1/0F); oxidizer flow rate through injector 
WT = Total propellant flow rate through injector 
p^x Oxidizer density at injector entry 
OF = Propellant mixture ratio 


242 


Since the computer program may be applied to a pump-fed engine with preburners, 
the equation for the pressure drop takes the following form: 

AP = WINJX^ (7_70) 

OX 

2.*g*Cp ^ DX^*NOX) ^ 

where 


WINJX = WOX - WXPREB 

(For definitions see Figures 2-10 and 2-11 in Section 2.4) . Finally, 

F32 = Ap (7-71) 

ox 

A flow chart for subroutines SECF31 and SECF32 and their interaction with 
Subroutine POSTAL is depicted in Figure 7-20. 


PRESSURE DROP^PARAMET 
TO SECFZO/ 21/ AND , 


r" 


■COMBUSTOR DESIGN! 
LOOP 


to 

o:oo 

UJCM 


POSTAL 


CALCULATES THE IN- 
JECTOR PRESSURE DROP 
FOR FUEL AND OXIDIZER1 
SIDE 


SECF31 


COMPUTES THE PRESSURE! 
DROP CHARACTERISTIC,. 

f31, for the fuel 

SIDE OF THE INJECTOR 


h 


T 


COMBUSTOR DESIGN 
LOOP 


SECF32 

COMPUTES THE PRESSUREl 
DROP CHARACTERISTIC, 
f32, for the OXIDIZER! 
SIDE OF THE INJECTOR 


CONTINUE COMBUSTOR 
DESIGN LOOP 


Figure 7-20. Flow Chart for Subroutines SECF31 and SECF32 

for Calculating the Pressure Drop Characteristics 
F31 and F32 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF3I AND SECF32 
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7.' 4 Combiisfor Complexity Characteristic' Modules (Subroutines)' 

7.4.1 Subroutine SECF41 . Subroutine SECF41 is the number of fuel plus 
oxidizer holes characteristic, and is given by 

F41 = (NF + NOX) - 1 (7-72) 

7.4.2 Subroutine SECF42 . Subroutine SECF42 is the volume of the 
oxidizer dome characteristic and is given by 

F42 = VOLX -■ (7-73) - 

7.4.3 Subroutine SECF43 . -Subroutine SECF43 is the volume of the fuel; 
dome characteristic and is given by 

F43’ = VOLF ■ (7-74) 

7.4.4 Subroutine SECF44 . Subroutine SECF44 is the length of the oxidizer 

holes characteristic and is given by : _ 

' . ’ >* • . 

F44-,= LOX • (7-75) 

7.4.5 S ubroutine SECF45 . Subroutine SECF45 is the length of the fuel 

holes characteristic and is given by ‘ -- , • 

F45 = LF (7-76)- 

7.4.6 Subroutine SECF46 . ^Subroutine SECF46 is the'* injector type com- 
plexity characteristic and is given by the index INJCMP (for definition, of 
INJCMP, see Section 2.4). At this time. Subroutine SECF46 is notlactive. 
Therefore , 

F46 = 0‘ . ’(7-77) 

However, INJCMP is an input parameter in Data Input Namelist COMMIS and 
is selected according to Figure 2-9, Section 2.4. 


Figure 7-21 depicts the flow chart for Subroutines SECF41 through SECF46. 
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Figure 7-21. Flow Chart for Subroutines SECF41 through 
SECF46 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF41 THROUGH SECF46 
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F46 -- I — 0.0 ‘ 0 Injector type complexity characteristic 



















• 7.5 Combustor Chamber Length Characteristic Module (Subroutine) 

I ■ . . ■ ' 

;• 7.5.1 Subroutine SECF51. Subroutine SECF51 computes the chamber length 

characteristic by 

I F51 = LC/DTHRTI (7-78) 

where 

I r LC = Combustion chamber length 

' DTHRTI = Diameter '.of throat 

Figure 7-22 depicts the flow chart for subroutine SECF51. 
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Figure 7-22. Flow Chart for Subroutine SECF51 
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Figure 7-23. iFlow Chart for Subroutine SECF61 
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.... .7.7 PropeHant. .Mixture Ratio Characteristic Module (Subroutine). 

! 

\ 7.7.1 Subroutine SECF71 . Subroutine SECF71 is the mixture ratii 

[teristic and is given by ^ 

j F71 = OF (7-80) 

I . 

^ OF as a free variable is perturbed by AESOP and the respective value 
.V in other subroutines. 

Figure 7-24 depicts the flow chart for Subroutine SECF71. 


COMBUSTOR DESIGN 
LOOP 


1 SEGF71 1 
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of' 
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CONTINUE COMBUSTOR 
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Figure 7-24. Flow Chart for Subroutine SECF71 


255 


charac 


used 
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■ 7.8 ‘ Baffle 'Chiafact'efistic Modules. (Subroutines) 


i 7.8.1 Subroutine SECF81. Subroutine SECF81 is a baffle characteristic 

! given by the number (integer) of radial vanes in baffle (BFLRAD) : 

1 F81 = BFLRAD (7-81) 

t 

? BFLRAD as an independent variable may be perturbed by Subroutine AESOP, 

^ and the respective value influences the specific combustor characteristics 
‘ F22, F23, F24, F25, F27, F28, and F29 as presented in Section 7.8.3. 

^ * 

’ 7.8.2 Subroutine SECF82. Subroutine SECF82 is, a baffle characteristic 

given by the number (integer) of circumferential cans in baffle (BFLCIR) : 

F82 = BFLCIR (7-?2) 

‘ BFLCIR as an independent variable may be perturbed by Subroutine AESOP, 
and the respective value influences the specific combustor characteristics 

• F22, F23, F24, F25, F27, F28, and F29 as presented in Section 7.8.3. 

1 . ' . ; * ' 

Figure 7-25 depicts the flow chart for Subroutines SECF81 and SECF82. 


7.8.3 Program Provisions to Accept Baffles (Subroutine FNDSNH) . The 
incorporation of baffles makes necessary the use of a new set of transverse 
■and radial frequency modes for Subroutines SECF24, 25, 27, and 28. The 
selection of these frequency modes is implemented by introducing proper 
exclusion rules. Subroutine FNDSNH serves this purpose. First, it contains 
the oscillation constants, s^}^, as functions of the transverse and radial 
mode numbers. Table 7-4 lists these values as they are contained in the 
.coded arrays. 

Second, it contains the exclusion rules for the selection of the lowest 
six s.^j^ values not forbidden by: 

V - F81/2, if F81 is ah even number 

V = F81, if F81 is an odd number, and 

h > F82, if F82 is an odd or even number . 

Note: For all s selections, the case where v= 0 with h > i, 

is alos included. 

Additional selection criteria for s.^^^ values are discussed in the respective 
sections for Subroutines SECF24, 25, 27, and 28. 
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Figure 7-25. Flow Chart for Subroutines SECF81 
and SECF82 
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PARAMETERS (COMMON) USED IN SUBROUTINE (S) SECF81 AND SECF82 
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TRANSVERSE MODE RADIAL MODE 
NUMBER- I number: h' 
(RADIAL BAFFLES) (CIRCUMFERENTIAL B.) 


0 

1 

1 

1 

2 

1 

0 

2 

3 

1 

4 

1 

1 

2 - 

- -5 

1 

2 

2 

0 

3 



S^l-VALUE 


0.00000 

1.8A118 

3.05A2A 

3.83170 

4.20h'9 


11.34592 
11.70600 
11.73494 
11.77088 
12.68191 > 


12.82649 

12.93239 

13.17037 

13.32369 

13.87884- 


TRANSVERSE MODE 
NUMBER. V 
(RADIAL BAFFLES) 




* # 

. 13 - 


7 


18 


10 


5 


■ 3 . . -■ - .. 


15.26818 
15.28674 
15 .-96411 




15.97544 

16.34752 

16.44785 

16.47063 

16.52937 


17.02032 
' 17.31284" 
'17.60027' 
17.77401 
17.78875 


18.01553 
- 18.06326 
18.63744 
18.74509 
19.00459 


19.10446 

19.19603 

19.51291 

19.61585 

19.88322 


19.94185 

20.14408 

20.22303 

20.57551 

20.97248 


21.01540 

21.16437 

21.18227 

21.22906 

21.43085 


21.93172 

22.14225 

22.21915 

22.40103 

22.50140 


22,62930 

22.67158 

22.76008 


■ 

I 

I 


TABLE; 7-4. Arrays for the Frequency Modes, s^j^, and their 
Components Stored in Subroutine FNDSNH 







































7.9 Minor Design Variable Dependent Functions 




There are several minor design variable dependent functions which are 
computed in Subroutine JPOSTAL (Section 6.6. 4) . Their analytical definitions 
^follow ..the ^iphabetical sequence ., i ,.-. - , .. - • . 

, . ■ ' ' ' 

7'. 9.1 AF, the Area of a Single 'Injector Fuel Orifice . 

J ; ? AF = 7I*DF**2/4: ' (7-83) 

where. DF is the, minimum diameter of the fuel orifice, and determines the 
pressure drop of the fuel side of the injectpr. (For more details see Section 
2.4, Figure 2-7,). ■ 


■’;.7.9.2 AI,- the Area of the Injector Face . 

AI = 7r*'bl**2/4 


(7-84) 


whefb.Dl = DC is the combustion chamber diameter at the, injector face. 

: 7*9.3 AOX, the Area of a Single, Injector Oxidizer Orifice . 

v''< „ ; ^ AOX = ■n*DX**2/4, ■ (7-85) , • 

where DX is the minimum diameter of the oxidizer orifice, and determines the 
pressure drop of the oxidizer side of the injectpr. (For more details see 
Section 2.' 4, Figure 2-7). 


7.9.4 CFIDEL, the; Ideal Thrust Coefficient . 

■ CFIDEL-=rTHRUST/(PC*ATHRTI) (7-86) 

where PC ■ combustion chamber pressure; ATHRTI = throat area. 

7.9.5 CS, the Speed of Sound in the Combustion Chamber . 

" CS' i Y fCOMB*SPi^T*RGAS*GREF ' (7-87) 


where 


TCOMB = Combustion temperature 
. SPHEAT = Speci’Hc heat ratio 
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... . RGAS =. Gas constant .per ..unit , weight 

GREF = Gravitational acceleration 

7,9.6 CSTR, the .Characteristic .Exhaust Velocity.' 


CSTR = 


CS 


SPHEAT\/,[2/CSPHEAT + 1]] 


' (SPHEAT ‘+ 1)/(SPHEAT -1) 


(7-88) 


where 


CS = Speed of .sound .in, combustion chamber 
SPHEAT = Specific heat ratio 


7.9.7 DELPF, the Pressure Drop Between Fuel Manifold and Combustion 
Chamber. 


where 


DELPF = 


,WINJF**2 


2*RHOF*GREF(CDF**2) ((NF*AF)**2) 

WINJF = Fuel flow rate for injector 
(see WINJF in this listing) 

RHOF = Density of fuel 

GREF = Gravitational acceleration 

CDF = Discharge coefficient for fuel orifices 

NF = Number of fuel holes 

AF = Area of a single fuel hole 
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7.9.8 DELPX-, the Pressure Drop’ Between Oxidizer Manifold and Combustion 
Chamber. 


where 


DELPX = 


WINJX**2 

'2*RH0X*GREF (CDX**2) ( (NOX*AOX>*’!'2}. . 


• (7-90).- 


WINJX = Oxidizer flow rate for injector 
C ee WINJX' In this listing) 

RHOX = Density of oxidizer 

GREF = Gravitational acceleration - 

CDX = Discharge 'coefficient for oxidizer orifices 
NOX = Number of oxidizer holes 


AOX = Area of a single oxidizer hole. 

' 

7.9.9 EPA, Injector Elements per Unit of Injector Face . 


. ,/EPA: =.NE/AI. ... . (7-91) 

where 


NE = Number of injector elements 
AI '= Area of injector face 


7.9.10 FJSPAC, Spacing of Fuel Holes . 


FJSPAC =■ ./ AI/NF 


(7-92) 


where 


Ai = Area of injector face 

NF = Number of fuel holes in injector 
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(7-93) 


7.9.11 ISP, Specific Impulse . For optimization runs-: 

ISP = YPOLYV(OF,PCISP,NPISP) 

For point design: 

PCISP(l) = — 

in Data Input Namelist CRVFIT (see Section 2. 7) . 

7.9.12 LD, Combustion Chamber Length/Diameter Ratio ; - ; 

LD = LC/DC (7-94) 

where 

LC = Length of combustion chamber 

DC = Diameter of combustion chamber at injector 
face 

LDF, Average Fuel Orifice Length/Diameter Ratio . 

LDF = LF/DF (7-95) 

LF = Length of the fuel hole 

(for definition see Section 2.3) 

DF = Diameter of the fuel hole 

(for definition see Section 2.3) 

7.9.14 LDOF, Log of Average Fuel Orifice Diameter ^- 

LDOF = ALOGIO(DOF) (7-96) 

' where 

DOF = 1000*DF 

DF = Diameter of the fuel hole 

(for definition see Section 2.3) 


7.9.13 


where 
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7.9.15 LDUM, Acoustical Length of Combustion Chamber. 


LDUM = LI + .667(L2 - LI) (7-97) 

where 

LI = Chamber length without converging section 

L2 = Chamber .including converging section (distance 
from injector face to throat area) 

(For more details see Section 7.2.3. 1). 


7.9.16 LRC, Log of Combustion Chamber Contraction Ratio. 


where 


LRC = ALOGIO(RC) 


(7-98) 


RC = Combustion chamber contraction ratio (for 
definition of RC see this listing) 


7.9.17 LSTR, the Characteristic Chamber Length . 

LSTR = VC/ATHRTI (7-99) 


where 


VC = Voliome of the combustion chamber (for 
definition of VC see this listing) 

ATHRTl = Throat area of the combustion chamber 


7.9.18 MC, Mach Number in Combustion Chamber 


MC = VELC/CS (7-100) 

where 


VELC = Average velocity of the gas in the combustion 
chamber (for definition of VELC see 
this listing) 

CS = Speed of sound in combustion gas 
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I 


7.9.19 
tion runs : 


/ 

( 


MOLWT, the Molecular Weight of the Combustion Gas. For optimiza- 


, MOLVfT = YPOLYV(OF,, PCMOLW, NPMOLW) 


For point design: 


PCMOLW(l) = 


in Data Input Namelist CRVFIT (see Section 2.7). 


(7-101) 


7.19.20 MPE, The Mass Flow Rate per Element . 

MPE = (WT. r, WFCP0L)/NE ,, (7-102) 


where 


WT = Total propellant mass flow rate 
WFCOOL =, Fuel flow rate for cooling 
NE = Number of injector elements 


7.9.21 NE, Number of Injector Elements . 

NE = (NF + NOX) /DIVIDE (7-103) 

where 

NF = Ni^ber of fuel orifices , 

NOX = Number of oxidizer orifices 
DIVIDE = 2 if INJCMP LE.3 
= 3 if INJCMP EQ. 4 
= 4 if INJCMP EQ. 5 
(values for INJCMP see Figure 2-9). 
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.7.9.22 PC, the Combustion, Chamber Pressure. 


PC = 


WT- 


ATHRTI*GREF*SPHEAT \ / [2/ (SPHEAT+1) ] 


(SPHEAT+ 1) / (SPHEAT- 1) 


GREF*SPHEAT*RGAS*TCOMB 


where 


(7-104) 


WT = Total propellant mass flows rate 
ATHRTI = Throat area- ^ 

- i . • ^ ' 

GREF = Gravitational constant 
SPHEAT = Propellant specific heat ratio 
RGAS = Propellant gas constant 
TCOMP = Combustion temperature 


7.9.23 PCOCRF, the Critical Pressure Ratio for Fuel. 


where 


PCOCRF^^ PCI/PCRITE .. ' (7-105) 

PCI = Chamber pressure 

PCRITF = Critical pressure of fuel (input for Input 
Data Namelist COMMIS) 


PCOCRF is used in Subroutines TAUNl and TAUN2 for computing the decay rates 
F,26 and F27. . - 
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7.9.24 PCOCRX, the Critical Pressure Ratio for Oxidizer. 



PCOCRX = PCI/PCRITX 

(7-106) 

where 

PCI = Chamber pressure 



PCRITX = Critical pressure of oxidizer (input for 
Input Data Namelist COMMIS) 

PCOCRX is used in Subroutines TAUNl arid TAUN2 for computing the decay rates 
F26 and F27. ' ■ 

7.9.25 

RC, the Combustion, Chamber Contraction Ratio. 



, RC = (D2/D3)**2 

(7-107) 

where 

D2 = Combustion chamber diameter 
(For a Cylindrical chamber: 

at nozzle entrance 
D2 = D1 =,DC) 


D3 = DTHRTI = throat diameter 


7i9.26 

RGASj the Gas Gpristarit per Unit Weight. 



RGAS = RGUNIV/MOLWr 

(7-108) 


where 

RGUNIV = Universal gas constant 

MOLOT = Propellant molecular weight; (for definition 
of MOLWT see this listing) 

7.9.27 SPHEAT, the, Specific Heat Ratio . For optimization runs: 
SRHEAT = YPOLYVCOF, PCK, NPK) (7-109) 
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For point design: ■ . . 

^ • '■ PCKCD = ■ 

in Data Input Namelist CRVFIT (see Section 2.7). 



7.9.29 TAX, the Total Area of Oxidizer Orifices. 

TAX = NOX*AOX (7-111) 

where 

• ’ NOX = ‘Number of oxidizer orifices 

AOX = Area of a single oxidizer orifice 



LSTR = Characteristic chamber length 

CSTR = Characteristic exhaust velocity 

i' - - ' ' 

SPHEAT = Specific heat ratio of the combustion gas 
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THETAG is used in Subroutines SECF20 and SECF21 (Section 7.2.1) for computing 
T E, J, and in Subroutine POSTAL (Section 7.9.34) for.,, computing VELC. Its 
definition is given in Reference 7-3, page 29. 

7.9.31. TCOMB, The Combustion Temperature. For optimization runs: 

■ TCOMB = YPOLYV (OF, PCTl, NPTl) " ' (7-113) 

For point design: /' ; ; ; 

. ■ ' PCTi(i) = 'v--- 

in Data Input Namelist CRVFIT (see Section 2.,7)- 


where 


7.9.32 TPVM, the Modified Thrust per Unit Volume . 

TPVM = PCI*ATHRTI/VC (7-114) 

‘ ■ PCI = Combustion chamber pressure 

ATHRTI^ = Throat area 

VC = Combustion chamber volume (for definition 
see this listing) 


7.9.33 VC, the Volume of the Combustion Chamber. 

VC = ~ (L1(Di2 + D2^ +D1*D2) + 

(L2 - L1)(D22 + D3^ + D2*D3)) 


(7-1 is) 
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where . 

LI * Length of the combustion chamber from the 
. „ injector, face to the onset of the converging 
section 

L2 = Length of the converging section plus LI 

bl = .Combustipri chamber diameter at injector face 

D2 = Combustion chamber diameter at the onset of the 
,, .. . converging, section. For cylindrical chamber 

it is Dl =• D2. • 

D3 = Throat diameter 


7.9.34 VELC, the Average Velocity of fhe Combustion Gases in the 
Combustion Chamber. 


VELC = LC/THETAG 


(7-116) 


where 


LC = Length of the combustion chamber from the 
injector face to the throat area 


THETAG = The mean gas ^ residence time (for definition 
see this listing) 


7.9.35 VFUEL, the Average Fuel Injection Velocity . 

VFUEL = WINJF/(NF*RH0F*Tr*DFVEL2/4) (7-117) 

where 




WINJF = Fuel mass flow rate passing through the 
injector ■ ' 

NF = Niimber of fuel orifices 

RHOF = Density of fuel (condition at injector . * 
face) 

DFVEL = Exit diameter of the fuel orifices 


7.9.36 VO, the Average Oxidizer Injection Velocity . 

VO = VOXID = WINJX/(N0X*RH0X*ir*DXVEL^/4) (7-118) 

where 

WINJX = Oxidizer mass flow rate passing; the' injector 
NOX = Number of oxidizer orifices 


RHOX = Density of oxidizer (condition at injector > 
exit) . 

DXVEL = Exit, diameter of the oxidizer orifices 


(7-119) 


' 7 -. 9 . 37 ' WF, the Total -Fuel Mass Flow .Rate ♦ 

WF = WT/1 + OF 

where ■ ' ■ - . r ' ■ v ■ I ^ . 

W =_Total propellant- mass flow rate 
OF = Propellant mixture ratio 

7.9.38 WINJF, the Fuel Mass Flow Rate Passing ' through the Injector. 


where 


•.WINJF = WF -t- WXPREB -/WFCOOL (7-120) 

■ i r -J .s'! .... 1 - • . - l " : . ' 

WF = Total fuel mass flow rate 

WXPREB = Oxidizer mass flow rate for preburner. If 
no preburner, WXPREB becomes 0 

WFCOOL = Fuel mass flow rate for cooling. If no cooling 
.V flow exists,.. WFCOOL becomes 0 


7.9.39 WINJX, The Oxidizer Mass Flow Rate Passing through the 
Oxidizer Orifices . ■ , . : 

. : . • WINJX = WOX - WXPREB (7-121) 


where . „ ,r . 

WOX = Total oxidizer mass flow rate (for definition 
• see this, listing) 

WXPREB = Oxidizer mass flow rate passing the preburner 
and therefore, injected with the fuel into the 
main chamber 
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